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ABSTRACT 



A radar target, acting as a scatterer of an incident 
electromagnetic wave, can be considered as a linear time- 
invariant system. Previous work has shown that the target's 
pole locations are independent of the incident electromagnetic 
excitation, including incident wave shape, aspect and 
polarization. This thesis develops the Kumaresan-Tuf ts and 
Cadzow-Solomon signal processing algorithms into computer 
routines and evaluates their pole extraction performance. 
Data used to evaluate the extraction algorithms includes 
synthetic and integral equation generated signals with 
additive noise, in addition to measurements of scattering by 
scale models made in an anechoic chamber. 
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I . INTRODUCTION 



A radar target, acting as a scatterer of a specified 
incident electromagnetic wave, can be considered as a single 
input, single output, linear time-invariant (LTD system for 
a fixed field observation point. The target can thus be 
considered as a transfer function with poles and zeros. Baum 
demonstrated at the Air Force Weapons Laboratory that a 
target's induced current response to an incident electro- 
magnetic wave has identifiable poles determined by the 
composition and structural geometry of the target [1] . In 
1974, Moffatt and Mains proposed that the target’s scattered 
field pole locations are independent of the incident 
electromagnetic excitation, including aspect and polarization 
[2]. Morgan has proven theoretically that, for the case of 
a conducting target, the scattering response contains complex 
natural resonances which are independent of the incident 
electromagnetic excitation [3] . By determining the poles of 
a target's response, aspect independent target identification 
can be accomplished through the use of electromagnetic natural 
resonances . 

Although the concept of radar target identification 
through the use of natural resonances was first proposed in 
1974 by Mains and Moffatt [2] , only recently have signal 
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processing techniques been applied to locate the poles in a 
radar target's response in the presence of noise. Kumaresan- 
Tufts [4] and Cadzow-Solomon [5] have each developed 
algorithms which have proven successful in the presence of 
noise. This thesis develops computer routines based upon 
these two algorithms and examines their respective performance 
and appropriateness using a variety of scattering data. 

A. THE PROBLEM 

Since the performance of signal processing methods varies 
under different conditions, a system employed to identify 
targets would possibly reach a decision based on the combined 
output of several signal processing methods. For example, the 
Kumaresan-Tuf ts and Cadzow-Solomon methods could be used to 
extract poles from the response of scale model targets. The 
information so gathered could be used to build a data base for 
comparison with data similarly obtained in actual field use. 
The results of this system would serve as one input to a 
larger system. Other methods would provide input to the 
system, such as the K-pulse method of Kennaugh [6] and the 
annihilation filter used by Dunavin [7] , Morgan and Dunavin 
[8] and Chen [9]. As the name suggests, an annihilation 
filter annihilates the target's poles. A system using the 
annihilation filter concept would contain many such filters, 
each previously designed to cancel the poles of a specific 



2 



known target. In actual field use, a radar target's response 
would be input into each of the filters, and the target 
selected would be that matching the filter whose output 
exhibits the lowest signal energy. 

A system used to identify radar targets would require the 
following concept of employment. First, information required 
by each of the sub-systems would be obtained for every target 
class of concern. In actual field use, this information would 
be compared against actual radar target responses. The system 
would then determine the identity of the target based on the 
input from each of its sub-systems. 

B . BACKGROUND 

Consider a perfectly conducting target illuminated by an 
electromagnetic field. The current induced on the surface of 
this target at a given point must satisfy the magnetic field 
integral equation (MFIE) , [10] 

j(r,0=2nxH'(r,0+JJ K(r,r ,0 d, 

^pv 

where n is an outward unit vector normal to the surface of 
the object, J is the surface current density, .H^. is the 
incident magnetic field, and K is a Green's function dyadic* 
The entire equation is most easily understood as the sum of 
driven currents and "feedback" currents corresponding to the 
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cross-product terir. and surface integral term respectively. 
The term driven by the magnetic f ield , 2nxH ^ forms the physical 
optics portion of the total current. Physical optics 
describes the cross-product term as the induced current 
without interaction with the rest of the body. The Green's 
function kernel describes the current at a point on the object 
due to the feedback of currents from every other point on the 
object, as previously illuminated by the incident field. The 
current at each point is then summed over the surface of the 
object. Note that the surface integral term is of principal- 
value type; the integral excludes the point f=f . 

Once the incident magnetic field is no longer present, 
the solutions of (1) are considered the natural modes of the 
object. These natural modes are of the form, J^exp(Sj^) . The 
natural resonance frequencies are of the form, 

Sn=°n+j^n (2) 

where 0^, is the damping rate in Nepers/sec and is the 
frequency in radians/sec. The natural resonances of (2) are 
functions of the structural geometry of the object and are 
independent of the incident magnetic field. To understand 
how these natural resonances are unique to the geometry and 
composition of the object, consider a set of points on the 
object previously illuminated by the incident field, so that 
H*=0. The current at a given point in the set is due to the 
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infinite number of feedback currents from every other point 
in the set. Recall that these feedbacks are described by the 
Green's function kernel in the integral term of (1). Since 
the set of points previously illuminated is physically located 
on the same object, the infinite number of paths that connect 
a point with all other points in the set is the same for all 
points in the set. The infinite number of paths are unique 
to the structural geometry of the object and correspond 
exactly to the infinite number of paths taken by currents 
which feedback to a given point via the Green's function 
kernel. Finally, the composition of the target determines the 
surface current density on the object. Although an infinite 
number of resonances exists in any object, only a limited 
number of these will be measurably excited by an incident 
field of finite bandwidth. These resonances described in (2) 
appear as complex conjugate pairs in the left-half portion of 
the s-plane. 

In the far-field, the back-scattered response of a target 
to an incident plane wave is of the form 



-s 






r-r 



/c)dS' 



( 3 ) 



where c is the speed of light and p is the unit vector whose 
direction matches that of the plane wave's propagation. 
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Equation (3) is the result of integrating the current at 
each point on the target surface for a fixed point in the far- 
field. Recall that the current at each point on the target 
is defined by (1). Thus, the back-scattered far-field can 
be obtained by substituting (1) into (3) : 



The currents in (1) produce the field in (4). In fact, each 
term in (4) corresponds to the term in (1) which produced it. 
Specifically, the first term in (4) describes the physical 
optics scattered field generated by the 2fixH* current which, 
of course, is the first term in (1). Similarly, the second 
term in (4) is produced by the source-free currents defined 
by the second term in (1) . Like the current described in (1) , 
the field in (4) is the sum of two terms, a driven term and 
a term containing feedbacks. 

The results of (4) can also be seen as two forms of the 
Singularity Expansion Method (SEM) developed by Baum [1] . As 
shown by Morgan [10] , during the early-time portion of the 
target’s response, the scattered field is composed of the 
physical optics scattered field and a "Class 2" form of the 
SEM expansion. The class 2 SEM expansion corresponds to the 
second term of (4) , wherein the coefficients are time- 
varying as the wave passes over the target, since the currents 




(4) 



6 



producing this portion of the field are integrated over a 
time-varying surface area. At the instant the wave passes the 
last point of the target, the physical optics field vanishes 
and the remaining term in (4) is produced by constant 
coefficients H . The coefficients H are constant at this 
instant since the surface area in the integral in (3) is now 
constant. This instant also marks the transition of (4) from 
a "class 2" SEM expansion of time-varying coefficients to a 
"class 1" SEM expansion of constant coefficients. The 
scattered field due to a plane wave is therefore composed of 
a physical optics term and a class 2 SEM expansion in the 
early-time, and a simple class 1 expansion in the late-time. 

Actual measurement of the scattered far-zone field would 
be greatly aided by knowledge of the transition time of the 
field from early time to late time. From [10], this 
transition for a monostatic radar would occur at At=T+2 (D+d)/c 
seconds after radar turn-on. Here, T is the pulse duration, 
D is the target's dimension along the direction of wave 
propagation, d is the distance between the target and the 
measurement point and c is the speed of light. 

The discussion presented in this section was extracted 
from work done by Morgan in [10] . The reader is referred to 
this work for a more detailed treatment of the material in 
this section. 
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C. HISTORY 

The results of the previous section form the basis for the 
hypothesis that the natural resonances found in the scattering 
response of a target to an incident electromagnetic wave are 
unique to that target. Additionally, only a finite set of 
these natural resonances are measurably excited by a wave of 
finite bandwidth. In 1974, Moffatt and Mains proposed that 
the extraction of resonances from a target's response to 
electromagnetic excitation could be used for target 
identification. This work related to earlier work in 1965, 
when Kennaugh and Moffatt first developed the concept of a 
radar target as a linear time invariant system. Poles in the 
z-plane are directly related to the natural resonances of a 
target 



where is given by 
seconds. Hence, 

identification. The 
discussed in the next 



(2) and At is the sampling interval in 
pole extraction involves resonance 
use of pole extraction algorithms is 
chapter . 
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II. POLE EXTRACTION ALGORITHMS 



The use of pole extraction algorithms to identify radar 
targets is discussed in this chapter. A brief discussion of 
two methods precedes the in-depth evaluation of the Kumaresan- 
Tufts and Cadzow-Solomon algorithms. The evaluation of the 
latter two algorithms occurs in two stages. First, each 
algorithm will be evaluated in its ability to extract poles 
from data with known poles. Some of the data processed was 
generated at various signal to noise ratios by a computer 
program written by Morgan [11] . Additional data was produced 
by Morgan's time-domain thin wire integral equation computer 
program [12]. In the second stage, a side by side comparison 
is made of poles extracted by each method using transient 
scattering measurements for a thin wire and for various model 
aircraft. Comparisons between the two methods are made as the 
aspect of the aircraft is varied. 

A. PREVIOUS WORK 

1. Direct Minimization 

The most direct way to determine the natural 
resonances in a target's response is to minimize the mean- 
square error between the modeled signal and the received 
signal. In [10], Morgan determined that the late-time target 
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response to a radar could be represented as a sum of damped 



sinusoids given by 



y(t) = 2 A 
1 = 1 



0,1 

e ' cos ((»),t+0,) 



( 6 ) 



The frequency, uj| , and damping rate, are the same 
parameters found in the natural resonance defined in (2). 
Phase, 6| , and amplitude. A,, are the remaining parameters. 
The representation in (6) is the sum of an infinite number of 
resonances. The sampled response to an incident wave of 
finite bandwidth can be modeled as 



y(n*t) 



N 



A,e 



0,nat 



COS (o),n6t+8, ) 



(7) 



where At is the sampling interval in seconds. The four 
parameters of (7) must be adjusted to minimize the sampled 
mean-square error signal 



®n=(y-yj 



( 8 ) 



between the actual discrete sampled received signal y and the 

n 

modeled signal . The processing required in this 
minimization problem is both inefficient and highly non- 
linear. Nevertheless, Chong used this method to process 
mathematically-generated data down to 15.0 dB signal-to-noise 
(SNR) ratio [13] . 
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2. Prony's Method 

As in direct minimization, Prony’s approach to 

resonance classification focuses on the late-time portion of 

a radar target's response. However, linear processing and 

root solving are used. The late-time response is modeled as 

the output of an LTI system of order Kj,.- Each signal received 

at some discrete sample, n, is considered to be the weighted 

sum of Kj, previous signals. Thus, the finite term 

approximation of the received late-time signal, y , is defined 

n 

by 






1=1 



n- 1 



(9) 



The z-transform of (9) is 



Kn ^ Kq-I . Kn-1 , 

z '’-b.z -b,z ° . . .-b. =0 



(10) 



'1“ *^2‘- Np 

The roots of this polynomial in z are the poles of the system 
model. Therefore, the key to extracting the poles in the 
system's response lies in solving for the coefficients b, of 
(9) . 

A set of Kj)+M received signals in M equations (9) can 
be arranged in matrix form as 



Kp-l 



Y • • • y 

■'m-I ■'K|3+M-2 



[!•] 



Yv 1 






( 11 ) 
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In Prony ' s original method, the data matrix is 
exactly determined, and the coefficient vector, b, is solved 
using linear computations. In the presence of noise, Prony 
overdetermines the data matrix by setting M>Kp and solves for 
the coefficient vector by obtaining the least-squares solution 
to the system of equations. 

The Prony method has two major problems. First, the 
poles obtained by the least squares solution to the 
overdetermined matrix may be strongly perturbed by noise [14] , 
since noise does not satisfy the causal model of the system. 
Second, the order of the system is generally not known a 
priori. When the estimated order is greater than the actual 
order, poles due to noise are generated. Prony ' s method 
offers no technique for distinguishing between the signal 
poles and the extra poles caused by overestimation of the 
system's order. If the estimated system order is less than 
the actual order, actual poles are lost and the remaining 
poles are perturbed from their true positions. 

B . KUMARESAN-TUFTS ALGORITHM 

The Kumaresan and Tufts pole extraction algorithm was 
developed by adapting Prony 's method to reduce the problems 
addressed in the preceding section. The Kumaresan-Tuf ts 
algorithm modifies the least-squares Prony method in three 
ways : 
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1. Processed signals are arranged in a data matrix 
based on a non-casual model of the system. 

2. The model of the system is deliberately 
overestimated , 

3. The system of equations determined by the above 
two criteria is solved by using singular value 
decomposition (SVD). 

Kumaresan demonstrates in [15] that the use of singular 
value decomposition tends to force the extra poles of the 
excess-order system inside the unit circle, while the non- 
causal arrangement of the signals tends to force the signal 
poles outside the unit circle. The excess order of the system 
model reduces the effects of noise on the actual poles. Since 
the noise is stationary and stable, it looks the same in 
forward and backward time. 

1. Equations 

Recall that in (9), Prony's technique defines the 
received late-time signal as the weighted sum of previous 
signals, where Kj, is presumed to be the order of the system. 
Kumaresan models the same late-time signal as the weighted sum 
of Kj, future signals, where is greater than the estimated 
order of the system. This non-casual model is given by 




( 12 ) 
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A system of M such prediction equations can be written in 
matrix form as 



r y, • • • y. n 






1 K[) 

• • 




SI 


• • 






_ yn- • -yvn-i _ 




L b'l J 



Or, in matrix notation, 




(13) 



As in Prony's method, the coefficients [b'l. are coefficients of 
a polynomial in z that models the system's late-time response. 
Two simple manipulations of either data matrix leads to the 
relationship between the coefficients of the Prony model and 
the prediction coefficients of the Kumaresan-Tuf ts model. 
With bQ=-i , a prediction coefficient is related to an 
autoregressive coefficient by 




(15) 



From the above relationship, it can be shown that the complex 
pole pairs of the causal model are merely conjugate 
reflections across the unit circle of the pole pairs in the 
non-causal model. 

2. Singular Value Decomposition 

The non-causal arrangement of late-time signals in a 
set of system equations, and subsequent processing through 
singular value decomposition, combine to separate the signal 
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and noise into orthogonal spaces. As discussed in the 
preceding paragraph, poles of the non-causal model are 
reflected outside the unit circle. Kumaresan demonstrates in 
[15] that the extra poles of the excess-order system can be 
forced inside the unit circle through the use of SVD. 

Singular value decomposition factors the i MXKp data 
matrix into the product of the matrices: 



The columns of U (MXM) are eigenvectors of 

j 

columns of V (KpXKp) are eigenvectors of ^y^y • If r is the 
rank of the data matrix, , the diagonal matrix z (MXKp) 
contains r singular values which are the square roots of the 
nonzero eigenvalues of both and Dy*^y.* rearranging 

the three matrices in the product, the pseudoinverse of ^y. can 
be obtained as 



Dy = Vi*U^ 



(17) 



where z"^ is a (KpXM) matrix whose singular values on the 
diagonal are the reciprocals of those in the I matrix. 
Finally, the coefficient vector b* of minimum Euclidian norm, 
is given by 




( 18 ) 
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The coefficient vector b* so obtained is the minimum length 
least-squares solution to (14). In other words, b* is the 
best possible solution to (14). In the case of noiseless 
data, the extraneous poles generated by the excess-order model 
will always be inside the unit circle when b*. is used. This 
result is generally true for noisy data. 

3. Bias Compensation 

Kumaresan and Tufts [4] observed that the addition of 
noise perturbed the singular values of the z matrix of (16). 
If the perturbation of these singular values is not 
compensated, both the signal poles and extraneous poles are 
biased towards the unit circle. Kumaresan and Tufts used a 
compensation method which reduced the bias in their work, but 
did not derive an analytical justification. In [16], Norton 
derived a more valid bias compensation method based on the 
eigenvalue shifting theorem. 

4. Kumaresan and Tufts Compensation 

If the actual order of the system is _, then the 
first Kp singular values of the Z matrix in (16) are non- 
zero. The remaining Kj,-Kp singular values are considered 
noise singular values and are zero in the case of noiseless 
data. The addition of noise perturbs the first Kp signal 
singular values and increases the noise to some non-zero 
value. Kumaresan and Tufts compensated for this increase in 
the singular values due to the noise by subtracting the 
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average of the noise singular values from the signal singular 
values. The noise singular values were then set to zero. 

5. Compensation Based on Eigenvalue Shifting Theorem 

As described in the previous section, the singular 
values of the matrix are the square roots of the 

eigenvalues of D„D^ and d^D • Assume the noisy data matrix 

y y- . y y ’ 

can be represented by Dy=S+Ni» where N is composed of the wide- 
sense stationary white noise process v,, given by 

(19) 

The expected value of DyD^ can be obtained by 




DyDj=E[ (S+N) (S+N)^]=E[sS^]+e[SN^]+E[NST]+E[NN^] (20) 



Since S is deterministic, E[SS^ ]=SS''^ . Assuming the noise is 
zero mean, the two cross products are zero. Because we assume 
the noise is wide-sense stationary and white, E[NN^]=o^I, 
where is the noise variance and I is the identity matrix. 
The expected value of 1 D thus becomes 

y y- 



E 




( 21 ) 
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similarly, the expected value of f the other source of 

singular values, is 

E[D^Dy]=S'^S+o2l (22) 

The assumption in the results of (21) and (22) is that the 
diagonals of E[N^N]=E[NN^] equals the noise variance .♦ 
Equations (21) and (22) show that in the mean, the squares of 
the singular values of are increased by the noise variance. 

The results lead to the method of eigenvalue 
compensation recommended by Norton in [16]. Recall from (16) 
that the eigenvalues of Dy are on the diagonal of the Z 
matrix returned by the singular value decomposition of Dy • 
If Kq' is the actual order of the system, and is the 

estimated order of the system then the remaining K^-Kp' 
singular values of the Z matrix can be squared and averaged 
to obtain an estimate of the noise variance, . These noise 
singular values can then be set to zero. The first 
singular values of the Z matrix are then squared and reduced 
by subtracting the estimate of the noise variance. The square 
root of the difference becomes the new first Kp singular 
values of the compensated Z matrix. Calculations according 
to (17) and (18) can then be carried out in a normal manner 
to obtain poles in the presence of the noise. Eigenvalue 
compensation requires an estimate of the actual order of the 
system. Methods to obtain this estimate are discussed in 
Chapter III. 
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6 . 



Performance 



The Kumaresan-Tuf ts algorithm was programmed in 
Fortran and tested on various types of data. The program 
appears in Appendix A. 

a. Synthetically Generated Data 

The starting point for evaluating the performance 
of the Kumaresan-Tuf ts algorithm was with synthetically 
generated data of the form given by (8) and shown here again 
for convenience 



N 



0,n6l 

e ‘ cos ((i),n6t+0, ) 



( 8 > 



Again, Aj , a, , o), , 0,, are the amplitude, damping rate, frequency 
and phase of a set of N damped sinusoids. Noisy data was 
created by adding stationary white noise. 

1. Noise Performance 

The algorithm was evaluated at various SNR's, 
ranging from 90.0 dB to 7.0 dB. These SNR's are ratios of 
signal energy to noise energy rather than the ratio of signal- 
to-noise power. Synthetic data so generated more closely 
resembles the exponential decay of signal power typical in 
actual radar measurements. 
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Figure 1 shows the signal produced by two s- 
plane poles at 90.0 dB. Figures 2 through 6 depict the poles 
extracted from this signal at SNR's ranging from 90.0 dB to 
7.0 dB . Obtained poles are shown at their positions within 
the upper right hand quadrant of the unit circle in the z- 
plane. Not shown are conjugates of each pole which are 
located below the real axis outside the figure boundaries. 

Figures 2 through 6 demonstrate outstanding 
performance on noisy data, even at SNR's of 7.0 dB . The 
scaling needed to show a discernible difference between 
results obtained at 30.0 dB and 7.0 dB would necessarily 
exclude one of the poles from the enlarged figure. The 
average distance of the trial poles obtained in the 7.0 dB 
SNR signal from the true poles is on the order of 10“^. This 
magnitude corresponds to that of the average estimate of the 
noise variance obtained in successive trials with this signal. 
The correlation between the distance of trial poles from true 
poles and the noise variance estimate was consistently 
observed with each of the different signal-to-noise ratios 
used. Figure 7 depicts the signal of Figure 1 severely 
corrupted by noise having 7.0 dB SNR. 

As discussed previously, the signal-to-noise 
ratio used in the synthetically generated data is the ratio 
of energy. Figure 8 depicts the results of pole extraction 
from the signal shown in Figure 7, but with a late-time 
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Damped Cosine: 4 s-plane poles; (-.l,+ -1.0),(-.4,+ -13.0) 90.0dB 
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Figure 1- 



Signal Containing two S-Plane Poles, 
90.0 dB SNR 
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512 samples over 20.0 nanosec 



Kumaresan-Tufts pole extraction 90.0dB 
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Figure 2 . Kumaresan-Tufts Poles, Synthetic Data, 90.0 dB SNR 
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Kumaresan-Tufts 30.0dB 
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Figure 3. Kumaresan-Tufts Poles, Synthetic Data, 30.0 dB SNR 
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Kumaresan-Tufts 20.0dB 




2 3eiuj 

Figure 4. Kumaresan-Tufts Poles, Synthetic Data, 20.0 dB SNR 
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Kumaresan-Tufts lO.OdB 
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Figure 5. Kumaresan-Tufts Poles, Synthetic Data, 10.0 dB SNR 
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Kumaresan-Tufts pole extraction 7.0dB 
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Figure 6. Kumaresan-Tufts Poles, Synthetic Data, 7 . 0 dB SNR 
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Damped Cosine: 4 s-plane poles; (-.1,+ -1.0), (-.4,+ -13.0) 7.0dB 
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Kumaresan-Tufts pole extraction, 7.0dB, early time 10.0 nanosecs 
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Figure 8. 



Kumaresan-Tufts Pole Extraction, 



7.0 dB SNR 
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Real z 



beginning ten nanoseconds later. Since the SNR is calculated 
over twenty nanoseconds for both signals, the signal power at 
some later time will clearly be less than the power ten 
nanoseconds earlier. The results in Figure 8 show complete 
breakdown of the algorithm's ability to extract poles. The 
trial poles shown are the poles closest to the true poles, and 
yet they are located at positions whose reflections are inside 
the unit circle where noise poles are typically located. 

The preceding results show outstanding 
accuracy for full-length noisy data but a complete breakdown 
of the algorithm for the same signal with a later transition 
to late-time. These initial observations are supported by 
similar findings presented in this thesis. 

b. Thin Wire Integral Equation Generated Data 

For simple objects such as a thin wire, the radar 
response of that object can be computed by establishing 
boundary conditions on the object and numerically solving the 
integral equations that describe the surface current. Recall 
the magnetic field integral equation given by (1). 
Simulations produced by Morgan's time-domain thin wire 
integral equation computer program [12] were used to evaluate 
the pole extraction algorithm. The excitation waveform used 
is the double Gaussian pulse depicted in Figure 9. This pulse 
is a wide Gaussian pulse with a ten percent width of 0.3 
nanoseconds subtracted from a narrow Gaussian pulse with a ten 
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Double Gaussian Curve 




Figure 9. Double Gaussian Pulse 
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nanoseconds subtracted from a narrow Gaussian pulse with a ten 
percent width of 0.15 nanoseconds. 

Figures 10 through 13 depict back scattering 
response of a 0.1 meter length thin-wire, having a radius of 
0.00118 meter, computed at various incident aspects, ranging 
from thirty degrees to ninety degrees. The laboratory 
arrangement for actual measurements simulated by Morgan's 
program is described in [17] . Ninety degrees represents a 
broadside aspect, while thirty degrees represents the incident 
plane wave having nearly grazing incidence on the wire. The 
poles extracted at each of the four aspect angles are plotted 
in Figure 14. In this figure, and those that follow which 
depict extracted poles, the signal poles lie in or on the unit 
circle, and the noise poles lie outside. 

The results obtained with this rigorous numerical 
computation demonstrate the aspect independence of the 
extracted poles using the Kumaresan-Tuf ts method. Note that 
only half of the poles were obtained for broadside 
illumination; two even-numbered poles can easily be seen 
outside the unit circle. This results because of the physical 
symmetry of both the wire and the incident field, thus 
precluding excitation of odd-symmetric modal currents and 
their associated natural resonances. 

Figure 15 exemplifies the computed back-scattering 
response of the 0.1 meter thin wire corrupted artificially 
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30 degree backscattering 
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Figure 10. Integral Equation Thin Wire Scattering, 30 Degree 
Aspect 
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45 degree backscattering 
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Figure 11. Integral Equation Thin Wire Scattering, 45 Degree 
Aspect 
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Figure 12. 



Integral Equation Thin Wire Scattering, 
Aspect 



60 Degree 
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90 degree backscattering 
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Figure 13. Integral Equation Thin Wire Scattering, 90 Degree 
Aspect 
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Figure 14, Kumaresan-Tufts Poles, Noiseless Thin Wire Data 
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Figure 15. 



Integral Equation Thin Wire Scattering, 20.0 dB 
SNR, 45 Degree Aspect 
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with noise at a 20.0 dB SNR. Figure 16 shows the poles 
extracted at each of the four angles of incidence used 
previously in Figure 14. Poles of Figure 14 at 90° are now 
missing in Figure 16, and only the first three low frequency 
poles are tightly grouped. The loss of high frequency poles 
is expected because these have the highest damping and thus 
lose their energy at the fastest rate. Further comparison 
between results computed at 20.0 dB SNR and infinite SNR are 
offered, angle by angle, in Figures 17 through 20. 

One additional test of the computed thin wire 
scattering was conducted at a 7.0 dB SNR. The corrupted 
waveforms are exemplified by Figure 21; the extracted poles 
are shown in Figure 22. The number of poles obtained has 
decreased with respect to the number obtained at 20.0 dB SNR. 
The grouping of the clusters has also expanded. Angle by 
angle comparisons are again offered in Figures 23 through 26. 
c. Scale Model Measurements 

The transient scattering measurements of scale 
models used for evaluation in this section were made by Walsh 
using the anechoic chamber of the Transient Electromagnetic 
Scattering Laboratory at the Naval Postgraduate School. The 
entire measurement process and laboratory setup are described 
in detail in [17] . 
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Figure 16. Kumaresan-Tufts Poles, 20.0 dB SNR 
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Figure 17 . 



Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 30 Degree Aspect 
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Figure 18. 



Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 45 Degree Aspect 
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Figure 19. Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 60 Degree Aspect 
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Figure 20. Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 90 Degree Aspect 
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Figure 21. Integral Equation Thin Wire Scattering, 7.0 dB 
SNR, 45 Degree Aspect 
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Figure 22. Kumaresan-Tufts Poles, 7 . 0 dB SNR 
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Figure 23. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 30 Degree Aspect 
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Figure 24. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 45 Degree Aspect 
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Figure 25. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 60 Degree Aspect 
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Figure 26. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7 . 0 dB SNR, 90 Degree Aspect 
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1 . 



wire Targets 



The thin wire measurements were obtained from 
the scattering response of a 0.1 meter length thin wire having 
radius 0.00118 meter. Recall that these are the same 
dimensions as the wire whose computed response was processed 
in the previous section. The measurements at each of four 
incident aspects are shown in Figures 27 through 30. 

The poles extracted from the four measurements 
are depicted in Figure 31. As before in the computed noisy 
data, tight clusters occur only at the lowest frequencies. 
The poles in these tight clusters are those which are 
measurably present at various aspects. The poles extracted 
at higher frequencies are those which possessed sufficient 
measurable energy at the given aspect. Figure 32 depicts the 
comparison between poles extracted from the measured and 
computed signals. Again, the closest agreement between the 
two sets of poles occurs at the lowest frequences. 

2. Aircraft Models 

Plastic 1/72 scale aircraft models, coated with 
silver, were used for transient scattering measurements. 
Representative scattering signatures of two aircraft targets, 
measured at six different aspects, are shown in Figures 33 
through 36. 

The results of pole extraction in target 1 are 
shown for a total of six different aspects in Figures 37 and 
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30 degree backscatiering 
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Figure 27. Measured Thin Wire Scattering, 30 Degree Aspect 
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45 degree backscattering 
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Figure 28. Measured Thin Wire Scattering, 45 Degree Aspect 
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Figure 29. Measured Thin Wire Scattering, 60 Degree Aspect 
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90 degree backscattering 
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Figure 30. 
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Measured Thin Wire Scattering, 90 Degree Aspect 
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Figure 31. 
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Kumaresan-Tufts Poles, 
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Figure 32. Thin Wire Comparison, Measured vs. 
Equation 
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Figure 33. Target 1 Scattering, 30 Degrees from Nose on 
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Figure 35. Target 2 Scattering, 30 Degrees from Nose on 



59 



Nose on 



xiO-3 




Figure 36. Target 2 Scattering, Nose on 
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Figure 37. Kumaresan-Tufts Poles, Target 1 
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38. The poles extracted at all six aspects are shown in 
Figure 39. Only one clearly discernible cluster is present 
in each of the three figures. At higher frequencies, no 
useful information is imparted by the data. Results of 
similar, though slightly improved quality, were obtained from 
target 2. These results are presented in Figures 40 through 
42 in the format of Figures 37 through 39 respectively. 

Although the Kumaresan-Tuf ts algorithm is 
capable of extracting low frequency poles acceptably, the 
inconsistent results at higher frequences reveals the inherent 
weakness in an algorithm capable of processing only the late- 
time portion of a target’s radar response. 

A side-by-side comparison of poles obtained from 
both aircraft by both the Kumaresan-Tuf ts method and the 
Cadzow-Solomon method is presented at the end of the chapter 
to illustrate the gains afforded by processing the early— time. 

C. CADZOW-SOLOMON ALGORITHM 

Recall from the results depicted in Figure 8 that a late 
transition to late-time, and the consequent reduction of 
signal power, caused complete breakdown of the Kumerasan-Tuf ts 
algorithm. The Cadzow-Solomon algorithm addresses this 
shortcoming by processing the signal at the instantaneous 
onset of early-time. Thus, the Cadzow-Solomon algorithm is 
capable of processing the earliest response of a target to 
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Figure 38. Kumaresan-Tufts Poles, Target 1 
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Figure 39. 



Kumaresan-Tufts Poles, Target 1 
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Figure 40. Kumaresan-Tufts Poles, Target 2 
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Figure 41. 



Kumaresan-Tufts Poles, Target 2 
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. Kumaresan-Tufts Poles, 



Figure 42 



Target 2 



electromagnetic excitation where the response has the greatest 
magnitude . 

1. Applicability 

The early-time portion of a target's scattered field 
occurs as long as there is a driven portion of the total 
field. Once the field no longer contains a scattered response 
due, in part, to the incident excitation at points on the 
object, early-time ceases and late— time begins. Hence, the 
Cadzow-Solomon models both the system's input and output, and 
equivalently, the poles and zeros of the system transfer 
function . 

2. Equations 

The Cadzow-Solomon algorithm extends the auto- 
regressive equation (9), used in Prony ’ s method, to the more 
general autoregressive moving average (ARMA) equation 

y=Ib,y +Sa,x„_, (23) 

n 1:, n-1 ,-o 

where the second summation term models the excitation to the 
system . 

A set of M such equations in matrix form is given by 



r y„ • • • y. . . . . x„ 



y • • *y Xw_i • • *Xi; +H-1 

Kp*M-2 ” ' 






y 



Kd-m-i 



(24) 
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As in the Kumaresan-Tuf ts method, M is selected to be greater 
than the column dimension of the data matrix which is 
Kd+Kn+1 . 

3. Excess Poles and Noise Removal 

The Cadzow-Solomon method used in this thesis is a 
modification which incorporates the non-causal arrangement of 
the system equations used by Kumaresan-Tuf ts . This 
modification was first discussed by Norton in [16] . The 
Kumaresan approach of overestimating the system order can be 
used as before in a non-causal model to constrain the noise 
poles inside the unit circle, while SVD forces the signal 
poles outside the unit circle. 

Since the input waveform is known, its order can be 
almost exactly determined. In all the work of this thesis, 
the input waveform used is the double Gaussian depicted in 
Figure 14. Approximately 25 samples defining this pulse of 
0.5 nanoseconds duration makes K^j equal 25 in equation (23). 
Since the input is causal, the signal zeros fall inside the 
unit circle where they cannot be easily segregated from 
similarly located noise poles. However, the signal zeros 
impart no information about the target and need not be 
extracted. The inclusion of the input in the data matrix is 
nevertheless vital to the model of the system and the accurate 
determination of the signal poles. 
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The ARMA equation of (23) can be modified to obtain 



U IX 

V = y b'.v +y^i>^r, I 

^ l=l ^ ‘ ‘ 1=0 



(25) 



The recursive portion of (25) is now in a non-causal form 
similar to expression (12) . A set of M such equations in 
matrix form is given by 












y • • *y +M-1 

l_ K^+M -^Kjj+Kq+M-1 ” * % ” ^ J 



b'v 



L “0 -1 



y. 



L _j 



Or, in matrix notation 



(26) 






(27) 



4. Singular Value Decomposition 

Lilce the system equations of the Kumaresan-Tuf ts 
model, the system equations in (26) are processed using 
singular value decomposition. The coefficient vector is again 
the minimum-norm solution, which constrains the extraneous 
poles and extraneous zeros to be inside the unit circle. 

5. Bias Compensation in the Cadzow-Solomon Formulation 
By compensating the eigenvalues of the i matrix in 

(16), the performance of the Kumaresan-Tuf ts algorithm is 
significantly improved in the presence of noise. Cadzow- 
Solomon have shown [5] that if the actual orders and 
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are 



overestimated to be and , min (Kp-Kp , K^-K^ ) singular 
values are zero in noiseless data. Since the input data is 
known, the eigenvalues of the data matrix may be compensated 
in the same manner as in the Kumaresan-Tuf ts algorithm for 
noiseless data. 

To understand the compensation required in noisy data, 
an analysis of additive noise is required. As given by Norton 
[16] , if the input data noise is w, and the output data noise 
is V, f the data matrix may be modeled as 
[Dy,] = [D :Dx]=S +N 

(28) 

where 

[Ny.] = [N,:N,] (29) 

and 
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The expected value of is then 

E [ ] =Sy,SJ,+E [ ] ( 31 ) 

If the input and output noise variances are not equal, the 
eigenvalue shifting theorem used in Kumaresan-Tuf ts cannot be 
used to analytically predict the requisite eigenvalue 
compensation of DyjjDyjj. Nevertheless, when the input and 

output variances were assumed equal, and eigenvalue 
compensation similar to that used in Kumaresan-Tuf ts was 
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performed, the results were consistently superior to those 
obtained without compensation. Therefore, the results of 
Cadzow-Solomon signal processing presented in this thesis were 
obtained using eigenvalue compensation and the assumption of 
equal noise variance. 

6. Performance 

The Cadzow-Solomon algorithm was programmed in Fortran 
and tested on the same data used for evaluating the Kumaresan- 
Tufts algorithm. Note that the Cadzow-Solomon algorithm 
can use the early-time portion of the data that the Kumaresan- 
Tufts algorithm can not use. The program appears in 
Appendix B. 

a. Synthetically Generated Data 

The starting point for evaluating the performance 
of the Cadzow-Solomon algorithm was with synthetically 
generated data of the form given by (8) plus the addition of 
input data required to model early time data. 

1. Noise Performance 

The algorithm was evaluated at various signal- 
to-noise ratios, ranging from 90.0 dB to 7.0 dB . Figure 43 
shows the signal produced by two s-plane poles at 90.0 dB , 
with a late-time beginning at 10.0 nanoseconds. Figures 44 
through 48 depict the poles extracted from this signal at the 
different signal-to-noise ratios. 
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Damped Cosine: 4 s-plane poles; (-.1,+ -1.0), (-.4,+ -13.0) SNR=90.0dB 




Figure 43. Signal Containing Two S-Plane Poles, 
90.0 dB SNR 



73 



512 samples over 20.0 nanosec Late time begins at 10.0 nanosec 



Cadzow-Solomon pole extraction 90.0db 
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Figure 44. Cadzow-Solomon Poles, Synthetic Data, 90.0 dB SNR 
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Cadzow-Solomon pole extraction 30.0db 
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Figure 45. Cadzow-Solomon Poles, Synthetic Data, 30.0 dB SNR 
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Cadzow-Solomon pole extraction 20.0db 
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Figure 46. 



Cadzow-Solomon Poles, Synthetic Data, 20.0 dB SNR 
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Cadzow-Solomon pole extraction lO.Odb 
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Figure 47. Cadzow-Solomon Poles, Synthetic Data, lO.O dB SNR 
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Cadzow-Solomon pole extraction 7.0db 
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Figure 48. 



Cadzow-Solomon Poles, Synthetic Data, 7.0 dB SNR 
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The figures chart the steady degradation of 
the algorithm's performance with the increase of noise. At 
30.0 dB, the location of the low frequency pole is already 
slightly displaced. More significant is the location of one 
of the extracted poles in the noise signal space. At 20.0 dB , 
the low frequency pole is located in some trials on the real 
axis. At 10.0 dB, all the extractions are located on the real 
axis and at 7.0 dB their locations there are dispersed. The 
extraction of the higher frequency pole is 
uncharacteristically more accurate than that of the low 
frequency pole. Even at 7.0 dB , the high frequency pole is 
located with excellent accuracy. The location of the low 
frequency pole near the real axis was chosen deliberately to 
illustrate the difficulty in resolving the slight frequency 
difference between the true pole and a noise pole located on 
the real axis. Also, fewer points were processed using the 
Cadzow-Solomon method than were processed using the Kumaresan- 
Tufts method, since the largest data matrix allowed by the 
programs in Appendices A and B contain fewer data points in 
the Cadzow-Solomon data matrix than in the Kumaresan-Tuf ts 
data matrix. The results demonstrate the need to process a 
substantial number of points in order to accurately extract 
low frequency poles. 
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b. Thin Wire Integral Equation Generated Data 

The performance of the Cadzow-Solomon algorithm 
was evaluated using the same set of data tested by the 
Kumaresan-Tuf ts algorithm. The results are presented in 
Figure 49. Tight clusters appear at frequencies higher than 
those obtained with the Kumaresan-Tuf ts algorithm. Figure 50 
depicts the poles extracted from the same signal at a 20.0 dB 
SNR. The clustering at this SNR is comparable to the results 
obtained by the Kumaresan-Tuf ts method with the noiseless 
data. Further angle-by-angle comparisons of the poles 
extracted from the noiseless data and the 20.0 dB data are 
depicted in Figures 51 through 54. Note the small number of 
poles in Figure 54 due to the unexcited odd-symmetric poles 
at 90° aspect. 

One further test was conducted on computed data 
at a 7.0 dB SNR. The results are depicted in Figure 55. Even 
at 7.0 dB , discernible clusters are present. Angle-by-angle 
comparisons of the poles obtained in 7.0 dB data and those 
obtained in noiseless data are presented in Figures 56 through 
59. 

c. Scale Models 

The same scale models used to evaluate the 
Kumaresan-Tuf ts algorithm were used to evaluate the Cadzow- 
Solomon algorithm. 
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Figure 49. Cadzow-Solomon Poles, Noiseless Thin Wire Data 
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Figure 50. Cadzow-Solomon Poles, 20.0 dB SNR 
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Figure 51. 



Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 30 Degree Aspect 
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Figure 52. Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 45 Degree Aspect 
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Figure 53. Integral Equation Thin Wire Comparison, Noiseless, 
vs. 20.0 dB SNR, 60 Degree Aspect 
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Figure 54 . 



Integral Equation Thin Wire Comparison, Noiseless 
vs. 20.0 dB SNR, 90 Degree Aspect 
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Figure 55. Cadzow-Solomon Poles, 7.0 dB SNR 
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Figure 56. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 30 Degree Aspect 
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Figure 57. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 45 Degree Aspect 
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Figure 58, 



tegral Equation Thin Wire Comparison, Noiseless 
7.0 dB SNR, 60 Degree Aspect 
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Figure 59 . 



Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 90 Degree Aspect 
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1. wire Targets 

Figure 60 depicts the poles extracted from 
measurements of a 0.1 meter wire. Three tight clusters appear 
at the lowest frequencies and at the highest frequencies. The 
poles in between can not be easily discriminated. The 
dispersion of these poles is apparently due to the aspect 
dependence of their measurable power. In other words, these 
poles are excited more at some aspects then at others. 

Figure 61 depicts the comparison between poles 
extracted from computed data and measured data. As in Figure 
60, close agreement exists at the highest and lowest 
frequencies. The results are much more favorable than those 
similarly obtained by the Kumaresan-Tuf ts algorithm. 

2. Model Aircraft 

Figures 62 through 64 depict poles extracted 
from aircraft target 1. As in the Kumaresan-Tuf ts testing, 
the Cadzow-Solomon testing was conducted at six different 
aspects. Results for target 2 are depicted in Figures 65 
through 67. The results of both targets show clearly defined 
clusters. The first two clusters of target 2 are 
exceptionally tight. However, the mid-frequency clusters of 
target 2 are not as clearly formed as those of target 1. 

Comparisons of poles obtained with each method 
for target 1 and 2 are depicted in Figure 68 and 69 
respectively. These two figures graphically depict the clear 
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Figure 60. Cadzow-Soloinon Poles, Measured Thin Wire 
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Figure 61. Thin Wire Comparison, Measured vs. Integral 
Equation 
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Figure 62. Cadzow-Solomon Poles Target 1, Three Aspects 
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Figure 63. Cadzow-Solomon Poles Target 1, Three Aspects 
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Figure 64. Cadzow-Solomon Poles Target 1, All Six Aspects 
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Figure 65. Cadzow-Solomon Poles Target 2, Three Aspects 
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Figure 66. Cadzow-Solomon Poles Target 2 , Three Aspects 
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Figure 67. Cadzow-Solomon Poles Target 2, All Six Targets 
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Figure 68. 



Pole Comparisons, Target 1, All Six Targets 
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Figure 69. 



Pole Comparisons, Target 2, All Six Aspects 
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superiority of the Cadzow-Solomon algorithm over the 
Kumaresan-Tuf ts algorithm. 

In order to obtain an initial indication of 
the possibility for target classification through pole 
extraction, nose-on measurements of two additional aircraft 
models were made, processed and compared with the results of 
targets 1 and 2. The nose-on measurements of targets 3 and 
4 appear in Figures 70 and 71 respectively. A comparison plot 
of poles extracted from each of the four targets is depicted 
in Figure 72. Each of the four aircraft measured are fighters 
of similar size and shape (see Table 1). The poles for each 
target are sufficiently different in this single measurement 
to identify each aircraft individually. However, some of the 
poles are arranged in clusters which appear with a harmonic 
pattern similar to that obtained for either of the first two 
aircraft at various aspects. In order to more fully assess 
the target classification capability of pole extraction, 
several measurements should be made of a given aircraft model. 
A plot of the poles extracted from each of these measurements 
would form clusters at the locations of the true poles. The 
centroid of each of these clusters would then be compared 
against the centroid poles similarly obtained from other 
aircraft. Although several poles of different aircraft might 
be similar, the set of poles belonging to an aircraft could 
form the basis for classification if that set was unique among 
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Figure 70. Target 3 Scattering, Nose-on 
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Figure 71. Target 4 Scattering, Nose~on 
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Figure 72. Cadzow- Solomon Pole Comparisons, 4 Targets, Nose-on 
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the sets belonging to all other measured aircraft. The results 
in Figure 72 demonstrate the possibility of using the Cadzow- 
Solomon pole extraction algorithm to aid in the classification 
of aircraft, perhaps by use of the extracted poles in 
constructing annihilation filters. 

TABLE 1. FULL SIZE DIMENSIONS OF TARGETS RECORDED 



Target number 


1 


2 


3 


4 


Overall length 
(meters ) 


12.20 


15.03 


16.94 


16.00 


Overall height 
(meters ) 


3.35 


5.09 


4.51 


4.80 


Wingspan 
(meters ) 


10.96 


10.00 


11 . 43 


13.95 


Tailplane span 


Unknown 


5.58 


6.92 


5.75 



(meters ) 
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III. SUMMARIES AND CONCLUSIONS 



In this chapter, a step-by-step guide through each 
algorithm is presented. At each step, techniques and lessons 
learned are discussed together with general observations. 
Conclusions are presented at the end of the chapter. 

A. KUMARESAN-TUFTS 

The first step in processing a signal with the Kumaresan- 
Tufts algorithm is to determine the beginning of early-time. 
The objective is to pick the earliest possible starting point 
without entering into the latter part of early-time. If the 
starting point for processing is improperly chosen to include 
the early-time, the results will be completely unreliable 
since the signal no longer satisfies the late time model. If 
the starting point is chosen too late, the signal may not be 
sufficiently strong in the presence of measurement noise. 
Since the signal is the sum of exponentially damped sinusoids, 
the optimum starting point is at the precise instant of 
transaction into late-time. The key to determining the 
beginning of late-time is in determining the beginning of 
early time. Determining the first response of the target to 
excitation cannot usually be done by a simple visual 
inspection of measurement data. Unless the exact distance to 
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the target is known, the most accurate method attempted by the 
author for determining the beginning of early— time is to 
process the signal using the Cadzow-Solomon algorithm. This 
is discussed in the next section. However, the reliance of 
the Kumaresan-Tuf ts algorithm on information provided by the 
Cadzow-Solomon algorithm is an obvious disadvantage of the 
former method. 

Once the starting point for processing has been selected, 
the next step is to determine the dimensions of the data 
matrix and, consequently, the number of points in the signal 
to be processed. In trials conducted on noiseless synthetic 
data, the accuracy of pole extraction increased steadily with 
the increase in the data matrix dimensions. These trials were 
conducted up to the limit of the array dimensions defined in 
the computer program of Appendix A. The number of points 
processed in measurement data should be as large as possible, 
while still meeting the following two constraints. First, 
incorporate as many cycles of the data as possible. Usually, 
visual inspection of the data reveals a repeating pattern 
which should be entirely incorporated into the window of 
points to be processed. When only portions of these patterns 
are selected, a disproportionate weighting tends to be placed 
on certain poles. Second, signal portions late in the 
response which are no longer distinguishable in the presence 
of measurement noise should not be selected. 
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The final step involves determining the number of true 
poles in the system. The following approach has proven to be 
the most successful. First, process the signal without any 
eigenvalue compensation to establish an upper bound on the 
order of the system. In most cases, the number of poles 
outside the unit circle will be less than the overestimated 
order of the system. If not, increase the row dimension of 
the data matrix in order to increase the estimated order of 
the system, and repeat. When the number of poles is less than 
the estimated order of the system, then one should gradually 
increase the number of eigenvalues compensated in successive 
trials, while closely observing the effects induced on the 
poles outside the unit circle. As the number of eigenvalues 
compensated is steadily increased, noise poles and weak signal 
poles will move inside the unit circle. The programs in 
Appendix A and B allow the user to compare the results of 
successive trials, by generating overlays for each plot. If 
N poles are in the signal space, at least the first N 
eigenvalues must not be compensated, or true poles may be 
lost. As the actual order of the system is approached by 
compensation, the user will notice an orderly, even 
arrangement assumed by the noise poles. If certain poles 
still remain suspect after compensation, vary slightly the 
other parameters, such as the starting point and the 
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dimensions of the data matrix. Generally, only true signal 
poles will repeatedly assert themselves under varying 
parameters . 

B. CADZOW-SOLOMON 

The techniques and general observations offered in the 
preceding section apply equally to the Cadzow-Solomon. 
algorithm. An important consideration in this method, not 
discussed above, is the selection of the beginning of early- 
time. Candidates for a starting point are usually at or near 
zero crossings within approximately thirty points of the 
object's first definite response to electromagnetic 
excitation. Begin processing at the chosen point while 
varying parameters in successive trials. Select the point 
whose successive results are the most consistent under varying 
parameters . 

The selection of the starting point for beginning of 
early-time can be very critical. For example, not a single 
pole could be extracted in one trial wherein the starting 
point occurred only ten points after the actual starting 
point. Additionally, in most cases observed, the late-time 
start given by the selected early-time occurred within less 
than two points from a zero crossing. If this observation 
proves to be generally true in later research, it may serve 
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as a way to check the starting point selected for one 
algorithm in terms of the other. 

C. CONCLUSIONS 

Both the Kumaresan-Tuf ts and the Cadzow-Solomon algorithms 
can effectively extract poles from the scattering response of 
a radar target. Because both algorithms obtain a least- 
squares solution to the system model, both perform acceptably 
in the presence of noise. Although eigenvalue compensation 
is not analytically justified in the Cadzow-Solomon algorithm, 
the results obtained through eigenvalue compensation in this 
method were generally superior to those similarly obtained in 
the Kumaresan-Tuf ts method. The results demonstrated the 
inherent advantages of an algorithm capable of processing a 
target’s strongest response in the early time. The Kumaresan- 
Tufts method compared favorably with the Cadzow-Solomon only 
in responses with a long late— time. 
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APPENDIX A. THE KUMARESAN-TUFTS POLE EXTRACTION ALGORTITHM 



The following program implements the Kumaresan-Tuf ts 

algorithm as described in Chapter 2 of this thesis. The 

program is written in Fortran 77. The SVD and root-finding 

subroutines called by this program are found in the EISPACK 

library [18] . The SVD subroutine is a translation from ALGOL 

as given in [19] . The matrix multiplication and graphics 

subroutines, also called by this program, are found in 

Appendix C and D respectively. 

INTEXIR im,Kd,M,l«,MAaai,NSTR1PT,I®LTRY 
lOTEm m,JOUS,NMENU,L/l/ 

IWIBGIR*2 KdPLT 

RERL*8 A(70,70) ,W(70) ,U(70,70) ,V(70,70) ,RV1{70) 

REAL*8 VS (70,70) ,UT(70,70) ,AINV(70,70) ,X(70) 

REAL*8 XP(70) ,B(70) ,SIGMA(70,70) ,SIG(70,70) 

REAL*8 CCF(70),R0arR(70),ROCrri(70) 

REAL*8 D(1024) ,AVG,MACHEP/1.0E-16/,Dy(140) 

C3CMPIiX*16 S(70) 

LOGICAL MATO/.TOUE./,MATV/.TOUE./,CAUSAL/.TOUE./,LDNG/.TRUE./ 

LOGICAL DSET/.FALSE./,NUFILE/.'TOJE./ 

CHARACTER TITIi^►16,HEADER*64,YN*l,IX:*l,T^ILER*16,TmZI*16 
CHARACTER TriL*16 

C Biter parameters for processing 

11 IF (DSET) CLOSE (10) 

NCVERLAY=0 

0PEN(10,FIli>='PLCT”) 

IF (DSET) GO TO 85 

HRTIE (*,*) 'Welccctie to signal processing using the' 

VRTIE (*,*) 'Kinnaresan-Tufts method' 

MOTE (*,*) ' ' 

MOTE (*,*) *Do you want ' 

MOTE (*,*) • ' 

MOTE (*,*) '1. The long version for beginners’ 

MOTE (*,*) '2. The short version for pros' 

MOTE (*,*) ' ’ 
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15 WRITE (*,*) 'Please enter 1 or 2 ’ 

READ (*,*) N 

IF (N .IQ. 1) 'IHEM 
LQNG=.TRUE. 

EliSEIF (N .IQ. 2) 1101 
I£W3=.FALSE. 

00 TO 15 
H©IF 

WRITE (*,*) 'Session vrill begin with entry of parameters needed fo+r processing' 
WRITE (*,*) 

WRITE (*,*) 'Do you want to enter parameters from' 

WRITE (*,*) ' ' 

WRITE (*,*) '1. The keyboard' 

WRITE (*,*) '2. A previously created file of parameters' 

WRITE (*,*) ' ’ 

16 WRITE (*,*) 'Please enter 1 or 2 ' 

READ (*,*) N 

IF (N .IQ. 1) TEm 
GO TO 1 

ELSEIF (N .IQ. 2) THEM 

10 WRITE (*,*) 'Biter title of file containirg parameters' 

READ (*,105) TTIL 
OPIN(l,FILE=Tnij) 

READ(1,105) TITLE 
READ (1,110) NPTS 
READ (1,110) NRT 
READ (1,110) Kd 
READ (1,110) M 
READ(1,110) DELTAY 
READ (1,110) NSTRTPT 
READ (1,110) NCAUS 
CLOSE (1) 

GO TO 85 
ELSE 
GO TO 16 
H©IF 

WRITE (*,*) ' ' 



1 NUFIL&.TRUE. 

IF (.NOT. DSET) NSTRTPT=1 

WRITE (*,*) 'Biter title of data file to be read' 

READ (*,105) TITLE 
OPBI(l,FIl£=Tni£) 

READ(1,105) HEAKR 
READ (1,110) NPTS 
IF (NPTS .CT. 1024) TWN 

WRITE (*,*) 'Nunber of points in data file exceeds the dimension' 
WRITE (*,*) 'of the array used in the program to store the file' 
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ST» 

n©iF 

aosE(i) 

IF (DSET) TOEN 

IF (NSIWPT+(Kd-HM)*DELTAy .!£. NPTS) GO TO 85 
H®IF 

3 IF (NUFIIi:) TOEN 

WRITE (*,*) 'Eiiter Kd, >= the estimated order of the system ' 
READ (*,*) Kd 
IF (Kd .CT. 69) TOEN 

WRITE {*,*) 'Kd must be less than 70, or dimension statements' 
WRITE (*,*) 'in this program must changed by the user' 

GO TO 3 

ELSEIF (Kd .LT. 2) TOEN 

WRITE (*,*) 'Kd must be at least 2' 

GO TO 3 
ENDIF 

IF (2*Kd .GT. NPTS) TOEN 

WRITC (*,*) 'Kd must be less than or equal to ',NPTS/2 
GO TO 3 

ELSEIF (2*Kd .EQ. NPTS) TOEN 
WRITO (*,*) 'Kd equals', Kd 
WRITE (*,*) 'M must be',Kd 
M=Kd 

WRITO (*,*) 'since there are a total of, NPTS 

WRITE (*,*) 'points in ',TnLE 

GO TO 45 

ENDIF 

GO TO 4 

ELSEIF (DSET) TOEN 

rwi 

20 IF (NSTRTPT+(N+M-1)*MLTAY .LE. NPTS) TOEN 

WRITE (*,*) 'Given the other parameters chosen thus far,' 

25 WRITE (*,*) 'Kd may range from ’,NRT 

WRITE (*,*) ' to',N 

WRITE (*,*) 'EiJter Kd' 

READ (*,*) Kd 

IF (Kd .(E. NRT .AND. Kd .LE. N) GO TO 85 
GO TO 25 
ELSE 
N^l 
GO TO 20 
ENDIF 
ENDIF 

4 IF (NUFIIE) TOEN 

WRITE (*,*) 'Ehter M, the row dimension of the data matrix' 

IF (.NOT. DSET .AND. LONG) TOEN 
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WRITE (*,*) ’ • 

WRITT (*,*) 'Note: Kd-m points in title 
WRITC (*,*) ' »dll be processed ' 

WRnr (*,*) • ' 

£}1DIF 

30 WRITE (*,*) *M may range fran',Kd 
IF (NPTSHa .(7T. 69) TffiN 
WRITC (*,*) ' to 69' 

ELSE 

WRITE {*,*) ' to',NPTS-Kd 

H®IF 

READ (*,*) M 
IF (M .Cr. 69) TOEN 

WRITE {*,*) 'M must also be less than 70' 

GO TO 30 

ELSEIF (M .LT. Kd) TOEN 

WRITE (*,*) 'M must be greater than or equal to Kd, Kd= ',Kd 
GO TO 30 

ELSEIF (Kd4M .err. NPTS) TOEN 

WRITE (*,*) 'Kd^M must be less than or equal to',NPTS, ',' 
WRITE (*,*) 'the number of data points in'.TITljE 
WRITE {*,*) ' ' 

GO TO 30 
ENDIF 
ELSE 
EH(d 

35 IF (NSrrRTPT+(Kd4N-l)*MLTAY .LE. NPTS) THEN 
EHtfl 
GO TO 35 
ELSE 
N=N-1 
ENDIF 

IF (N .EQ. Kd) THEN 

WRITE (*,*) 'M must equal', Kd 

M=Kd 

GO TO 85 

ENDIF 

IF (N .<?r. 69) 1^9 

40 WRITE {*,*) 'M may range from', Kd 
WRITE (*,*) ' to',N 

WRITE (*,*) 'Eiiter M' 

READ (*,*) M 

IF (M .GE. Kd .AND. M .I£. N) GO TO 85 

GO TO 40 

ENDIF 

45 IF (.NOT. NUFII£) GO TO 85 

5 Et=l 

50 IF (NSTRTPT4N*(Kd-HM) .LE. NPTS) THEN 
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55 



00 TO 50 
ELSE 
N=4f-1 
EM)IF 

IF (N .BQ. 1) TOEN 

WRITE (*,*) 'Given the other parameters chosai thus far,' 
WRITE (*,*) 'Spacing can only be 1' 

DELTRY=1 

IF (NUEILE) THEN 
GO TO 60 
ELSE 
GO TO 85 
ENDIF 
ENDIF 



IF (.NOT. DSET .AND. LONG) THEN 
WRITE {*,*) 'Ehter spacing between the ',Kd-41 
'data points of ' .TTHE 
'to be processed ' 



WRITE (*,*) 

WRITE (*,*) 

WRITE (*,*) 

WRITE (*,*) 

WRITE (*,*) 

WRITE (*,*) 

WRITE (*,*) 

ENDIF 

WRITE (*,*) 

WRITE (*,*) ' 

READ (*,*) MLTAY 
IF (DELTAY .GE. 1 
IF (NUFI1£) THEN 
GO TO 60 
ELSE 
GO TO 85 
ENDIF 
ELSE 
GO TO 55 
ENDIF 



'If, for example, one is chosen, then ',Kd+M 
'consecutive points in ',Tmi: 

'will be processed ' 



'Spacing may range fran 
' to',N 



.AND. DELTAY .LE. N) THEN 






60 WRITE (*,*) 'Do you wish to adjust eigenvalues? (y/n) ' 

READ (*,120) YN 

IF (YN .IQ. 'N' .OR. YN .IQ. 'n') THEN 
IF (NUEILE) GO TO 6 
GO TO 85 
ENDIF 

IF (YN .NE. 'Y' .AND. YN .NE. 'y') GO TO 60 
2 WRITE (*,*) 'Discard or ccnpensate eigenvalues? (d/c) ' 

READ (*,120) DC 

IF (DC .IQ. 'D' .<*. DC .IQ. 'd') GO TO 65 
IF (DC .NE. 'C .AND. DC .NE. 'c') GO TO 2 
WRITE (*,*) 'Eiiter estimate of the actual order of the system' 
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vRm (*,*) ’ ' 

IF (LCNG) THOI 

WRITE (*,*) "Riis estinate will be used to determine the ' 
WRITE (*,*) 'number of eigenvalues con^iensated or discarded ' 
ENDIF 

65 WRITE (*,*) 'the estimate may range fran 2' 

WRITE (*,*) ' to',Kd-l 

READ (*,*) NRT 

IF (NRT .or. Kd .CR. NRT .LT. 2) THEN 
00 TD 65 

ELSEIF {.NOT. NUFIIZ) THEJI 

GO TO 85 

ENDIF 



6 N?TRTTT=1 

70 IF (NSTOITT+(Rd-fM-l)*raLTAY .!£. NPTS) THEN 

NSTRTTT=f(STRTPT+l 
GO TO 70 
ELSE 

KSTRTPW«STR7TT-1 

ENDIF 

IF (NSTRITT .BQ. 1) TNEN 

WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE {*,*) 'the starting point for processing the data' 
WRITE (*,*) 'must be the first point in the data file' 

00 TO 85 
ENDIF 

WRITE {*,*) 'Ehter desired starting point in data file' 

IF (.NOT. DSET .AND. LONG) THEN 

WRITE {*,*) '1 indicates the first point in the data file ' 
ENDIF 

WRITE (*,*) ' ' 

WRITE (*,*) 'Givai the other parameters chosen thus far, ' 

75 WRITE (*,*) 'the starting point may range frcm 1' 

WRITE (*,*) ' to',NSTRTPT 

READ (*,*) N 

IF (N .GE. 1 .AND. N .!£. NSTRTPT) TEEN 

NSTETPT=N 

ELSE 

WRITE {*,*) 'Ehter starting point again' 

WRITE (*,*) ' ' 

GO TO 75 
ENDIF 

IF (.NOT. NUFILE) 00 TO 85 



7 WRITE (*,*) 'Do yew want the data matrix arrangemait to be' 
WRITE (*,*) ' ' 

WRITE (*,*) '1. Causal' 

WRITE (*,*) '2. Ncxi-causal' 

WRITE (*,*) ' ' 
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80 WRITE (*,*) 'Please enter 1 or 2 ' 

READ (*,*) NCAUS 
IF (NCAUS .BQ. 1) TOEK 
CAUSALf.TOUE. 

tLSEIF (NCAUS .EQ. 2) TOEN 
CAUSAL=.FALSE. 

ELSE 
00 •TO 80 
EIIDIF 
GO TO 85 

9 WRITE (*,*) 'Eiiter title of file to contain parameters' 
READ (*,105) TITL 
CfOl(l,FTLE!==TTTL) 

WRITE (1,105) TITliE 
WRITE (1,110) NPTS 
WRITE (1,110) NRT 
WRITE (1,110) Kd 
WRITE (1,110) M 
WRnE(l,110) MLTAY 
WRITE (1,110) NSTETPT 
WRITE (1,110) NCAUS 
CLOSE (1) 

IF (DSET) GO TO 85 

12 IF (DSET) TEEN 
CL0SE(2) 

CLOSEO) 

CAUi SUBPLT(NOVENLAY) 

ENDIF 
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DSin^.TEUE. 



NUFILE^. FALSE. 



WRITE(*,*) 

WRITE(*,*) 

+m£ 

WRITE(*,*) 

WRITE(*,*) 

WRITE(*,*) 

WRITE(*,*) 

WRITE(*,*) 

+Y 



'1. Data file to be processed ',T 

' Number of data points in data file ',NPTS 

'2. Estimated order of the system ',NRT 

'3. Kd, the number of colunns in the data matrix', Kd 
'4. M, the number of rows in the data matrix', M 
'5. Spacing between data points being processed ',MLTA 



WRITE(*,*) '6. First point in the data file to be processed' ,NSTET 
+PT 

WRITE(*,*) ' Last point in the data file to be processed' ,NSTET 



+PT+KiWM 

IF (NCAUS .BQ. 1) TEEN 

WRITE (*,*) '7. Data matrix arrangement for processing CA 

•KJSAL ' 

ELSE 



IF 
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NCK-CA 
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WRITE (*,*) '7. Data matrix arraiigement for processing 
-KJSMi ’ 

ENDIF 

WRITE(*,*) ’ ' 

WRITE(*,*) '8. Begin processing using above settings* 

WRTn:(*,*) *9. Store parameters 1-7 in a file* 

WRrra:(*,*) *10. Retrieve parameters 1-7 from a previously created 
+file* 



WRTITK*,*) *11. Reset overlays' 

WRrn:(*,*) *12. Re-plot overly* 

WRiraU*,*) *13. Bid this session of Rumaresan-Tufts signal process 
+ing* 

WRITE(*,*) ' ' 

WRrra:(*,*) 'Enter an integer from 1 to 12 to maXe dianges as often 
+ as you desire' 

READ (*,*) NMEIW 

IF (NMENU .LT. 1 .CR. NMEWU .CT. 13) 7®!! 



WRTIE(*,*) 'Biter an integer from 1 to 13' 

00 TO 90 

D©IF 



GO TO (1,2, 3, 4, 5,6, 7, 8, 9, 10, 11, 12, 13), NMENU 

8 CPIN(l,FII£==nTLE) 

R£AD(1,105) HEADER 
READ (1,110) NPTS 
READ (1,115) XQ 
READ(1,115) XQ 
DO 95 I=1,NPTS 
READ(1,115) D(I) 

95 CCWITNUE 

CLOSE (1) 

KdPL'WOd 

WRITE (*,*) 'Biter title of file to contain real part of poles' 

READ(*,105) Tnm 

0PiN(2,file=nTUR) 

WRITE(*,*) 'Biter title of file to contain imaginary part of poles' 

READ(*,105) TITLEI 

CPEN(3,file=rnLEI) 

WRriE(10,100) (KdPLT) 

WRITE(10,105) TnUF 
WRITE (10, 105) TTILEI 
100 FORMAT (12) 

MTH1AX(M,Rd) 

105 FORMAT (A) 

110 FC»MAT(I5) 

115 F»MAT(E12.6) 
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FORMAT (Al) 



C Form data matrix 

DO 125 I=l,Kd4M 

Dy (I)=D( (I-l) *DJLTAY4NSTOTT>T) 

125 OCKTINUE 

130 DO 140 1=1, M 
DO 135 J=l,Kd 
A(I,J)=Dy(I+J) 

135 OCKTIMJE 
140 CONTINUE 

B(l)=Dy{l) 

DO 145 1=2, M 
B(I)=A(I-1,1) 

145 CONTINUE 

C Begin singular value decanposition 

CALL SVD(MACHEP,M,Kd,MN,A,W,MA1U,U,MATV,V,im,RVl) 

C Errors in SVD? 

IF (HRR .(?T. 0.0) 'IHEN 

WRITE (*,*) 'Error in singular value number ',n3®,STC*> 
n©iF 

IF (YN .IQ. 'N') GO TO 190 

DO 150 1=1, Kd 
XP(I)=0.0 
150 CONTINUE 

C Discard or ccnjjensate eigenvalues 

C Order singular values 

XP(1)=W(1) 

DO 165 1=2, Kd 
DO 160 J=1,I 

IF (W(I) .<7T. XP(J)) TOEN 
DO 155 K=I+1,J,-1 
155 XP(K)=XP(K-1) 

XP(J)=W(I) 

GO TO 165 
EWDIF 

160 CONTINUE 

XP(I+1)=W(I) 

165 CONTINUE 

C XP{ ) now contains ordered singular values-XP(l) is the largest 
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C Discard eigenvalues 

IF (DC .BQ. 'D') nffiU 
DO 170 J=«RT+l,Kd 
170 W(J) = (0.0) 

ELSE 

C Conpensate eigenvalues 

kVO=O.Q 

DO 175 J=NRTfl,Kd 
AVG=AVOfXP{J)**2 
175 OCNTINUE 

IF (Kd .<?T. NRT) AVG=AVG/I»LE(FI£AT(KdH«T)) 

DO 185 J=l,Kd 
DO 180 K=l,Kd 

IF ( W(J) .BQ. XP(K) ) TWEN 
IF ( K .CT. NRT ) TOEN 
W(J)=0.0 
ELSE 

W(J)=DSg?T(DABS( W(J)*W(J)-AVG)) 

EJ©IF 
GO TO 185 
ENDIF 
180 CCKTINUE 
185 OOWnNUE 
EM)IF 

190 DO 200 1=1, M 
DO 195 J=1,M 
UT(I,J)=(U{J,I)) 

195 OCNTINUE 
200 OCWTINUE 

c Form SIGMA+ (KdxM) 

DO 210 1=1, Kd 
DO 205 J=1,M 
SIGMA(I,J)=0.0 

IF (I .EQ. J .AND. W(J) ,NE. 0.0) THEN 
SIGMA(I,J)=1.0D0/W(J) 

ELSE 

SIGMA(I,J)=O.OdO 

ENDIF 

205 OCNTINUE 
210 OCNimiE 

C Form SIGMA (Mrf(d) 

DO 220 1=1, M 
DO 215 J=l,Kd 
SIG(I,J)=0.0 

IF (I .BQ. J) SIG{I,J)=4f(J) 



122 



215 CXWTINUE 
220 OCNTIMJE 

C V=^&3rfOi,SI(>lA+^GiSxM,VS=4Ci3^ 

OIL !CMUL(V,SIGMA,Rd,Kd,M,VS) 

C VS=KdxM,Ur=MxM,AlNV=«d^ 

OUi MXMUL(VS,m’,Kd,M,M,AINV) 

C Calculate matrix multiplication of ADJV x B, where 

C AINV=KdxM,&^,XP=«dxl 

CAIL MXMUL(AINV,B,Kd,M,L,XP) 

C Calculate autoregressive coefficients from prediction coefficients 

IF (XP(Kd) .EQ. 0.0) raiN 
WRITE (*,*) avoiding division by zero' 

ELSE 

B(Rd)=1.0d0/XP(Kd) 

EI®IF 

DO 225 1=2, Kd 
B (I-l) =-« (Kd) *XP (Kd-I+1) 

225 CCNITNUE 

DO 230 1=1, Kd 
X(I)=-B (Kd-I+1) 

IF (NCAUS .EQ. 1) X (I) =-XP (Kd-I+1) 

230 ocwnuuE 
X(Kd+l)=1.0 

C Compute the roots of the polynomial in z 

CAUi P<MT(X,CCF,KD,ROCTO,ROan,IIR) 

IF (m .NE. 0) WRITE (*,*) 'E3«C» with PCXJ^T, m=’ ,IE2(,ST»> 

DO 235 1=1, Kd 
WRITE(2,115) ROCTO(I) 

WRITE(3,115) ROCTI(I) 

S (I) =DCMPLX (ROCTO (I) ,ROCrn (I) ) 

235 OOWTINUE 

MAGPCarO 
DO 240 1=1, Kd 

IF (CDABS(Sd)) .GE. l.OdO) MA0P(Xr=MA(3a/+l 
240 CCNTINUE 

WRITE (*,*) '# of poles with magnitude <= l',Kd-MR(3KXj 
WRITE (*,*) 'HIT ANY KEY TO CCNTINUE' 

READ (*,105) HEADER 
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C Plot poles 

NCraiIAY=4«JVHUAY+l 

CL0SE(2) 

CLOSEO) 

CaUi SUBPLT(NCVEFIAY) 

J=0 

K=0 



DO 245 1=1, Kd 

IF (CDABS(S(I)) .LT. 1.0) TOEN 

J=J+1 

K=K+1 

WRITE (*,*) S(I),0)ABS(S(I)) 

BNDIF 

IF (J .B3. 20) TOEH 

WRITE (*,*) 'Eiiter any key to continue' 

READ (*,105) HEM®3( 

H®IF 

245 CCMTINUE 

WRITE(*,*) 'Poles with magnitude less than one: *,K 

GO TO 85 
13 STCF 
M) 
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APPENDIX B: THE CADZOW-SOLOMON POLE EXTRACTION ALGORITHM 



The following program implements the Cadzow-Solomon 
algorithm as described in Chapter 2 of this thesis. The 
program is written in Fortran 77. The SVD and root-finding 
subroutines called by this program are found in the EISPACK 
library [18] . The SVD subroutine is a translation from ALGOL 
as given in [19] . The matrix multiplication and graphics 
subroutines, also called by this program, are found in 
Appendix C and D respectively. 



Slarge; 

INIBOTR nRR,M,Kn,M,m,MAGP(X,NSTRTPT,I®LTAY 
INTEGER m,NCAUS,NMaW,INSTRTPT 
INIBCaR*2 KdPLT 

REAL*8 A(70,70) ,W(70) ,U(70,70) ,V(70,70) ,RV1(70) 

REAL*8 VS(70,70) ,irT<70,70) ,AINV(70,70) ,X(70) 

REAL*8 XP(70) ,B(70) ,SI01A(70,70) ,SIG(70,70) 

REAL*8 OCf'(70),R0aiR(70),ROOTI(70) 

REAL MAG 

REAL*8 D(1024) , AVG,MACHEP/1.0E-16/, 1^(140) ,Dx(1024) 

0CMPLEX*16 S(70) 

IDGICAL MATU/ .TRUE. / ,MATV/ .TRUE. / ,CAUSAL/ .TRUE. /, LONG/ .TRUE. / 
LOGICAL DSET/.FALSE./,NUFILE/.TRUE./ 

CHARACUR TTILE*16,HEAI)HR*64,YN*l,I)C*l,TTTLER*16,TTim*16 
CHARACTER TTrL*16,TTni>*16 

C Biter parameters for processing 

14 IF (DSET) CLOSE (10) 

NCMRLAY=0 

cpEH(io,nLB='PLar’) 

IF (DSET) GO TO 215 

WRITE {*,*) 'Welccme to signal processing using the' 

WRITE (*,*) 'Cadzow-Solcnon method' 

WRITE (*,*) ' ' 

WRITE (*,*) 'Do you vrant ' 
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WRITE (*,*) ’ ’ 

WRITE (*,*) '1. The Icxig versicm for beginners' 
WRITE (*,*) '2. Tbe short version for pros' 
WRITE (*,*) ' ' 

25 WRITE (*,*) 'Please enter 1 or 2 ' 

READ (*,*) N 
IF (N .EQ. 1) THEN 
I£WG=.TOUE. 

ELSEIF (N .IQ. 2) TBEW 
LONG=.FALSE. 

ELSE 
GO TO 25 
EMDIF 



WRITE (*,*) 'Session will begin with entry of parameters needed fo 
+r processing' 

WRITE (*,*) 

WRITE (*,*) 'Do you want to enter parameters from' 

WRITE (*,*) ' ' 

WRITE (*,*) '1. The keyboard' 

WRITE (*,*) '2. A previously created file of parameters' 

WRITE (*,*) ' ' 

35 WRITE {*,*) 'Please enter 1 or 2 ' 

READ (*,*) N 
IF (N .EQ. 1) IHEN 
GO TO 8 

ELSEZF (N .EQ. 2) TBEJI 

13 WRITE (*,*) 'Eiiter title of file containing parameters' 

READ (*,100) TUL 
CPEN(l,FILE=nrL) 

READ(1,100) TIHE 
READ (1,110) NPTS 
READ (1,110) NRT 
READ (1,110) Kd 
READ(1,110) M 
READ(1,110) DELTAY 
READ (1,110) NSTRTPT 
READ (1,110) fOUS 
READ(1,100) TTriiD 
READ (1,110) MOTS 
READ(1,110) Kn 
READ (1,110) INSTRTPT 
CLOSE(l) 

GO TO 215 
ELSE 
GO TO 35 
ENDIF 

WRITE (*,*) ' ' 



8 WRITE (*,*) 'ELter title of file containing excitation waveform' 
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READ (*,100) TTITD 
CPEN(8,FILE=TnU)) 

READ(8,100) HEADER 

READ(8,U0) N 

IF (N .CT. 1024) TBEN 

WRITE (*,*) 'Number of points in data file e?cceeds the dimension' 

WRITE (*,*) 'of the array used in the program to store the file' 

STOT 

ENDIF 

CU0SE{8) 

IF ((N .GE. NDPTS) .AND. DSET) THEN 

NDPTS=ei 

GO •TO 215 

H®IF 

NIFrS=N 

9 WRITE (*,*) 'Enter estimated order of waveform' 

IF (DSET) THEN 

MAXIMLIM=NDPTS^! 

IF (MAXIMUM .GT. M-Kd-1) MAXIMUIf=4H(d-l 

IF (MAXIMUM .GT. NDPTS-INSTRTET-Kn-M+1) TEEN 

MAXIMUM^PTS-INSTOTFT-Kn-Mfl 

ENDIF 

ELSE 

MAXCMU»=66 

ENDIF 

IF (MAXIMUM .BQ. 1) TEEN 

WRITE (*,*) 'The estimated order of the waveform can only be 1' 
IF (DSET) GO TO 215 
GO TO 10 
ELSE 

IF (DSET) TEEN 

WRITE (*,*) 'Given the other parameters chosen thus far,' 

ENDIF 

45 WRITE (*,*) 'the order may range from 1' 

WRITE (*,*) ' to', MAXIMUM 

READ (*,*) Kn 

IF (Kn .($. 1 .AND. Kn .I£. MAXIMUM) TEEN 
IF (DSET) GO TO 215 
GO TO 10 
ENDIF 

WRITE (*,*) 'Eiiter estimated order again' 

WRITE (*,*) ' ' 

GO TO 45 
ENDIF 

IF (DSET) GO TO 215 

10 INSTETTT=1 

55 IF (INSTETPT+Kn4M-l .CT. NDPTS) TEEN 
INSTRTTT=INSTRT?T-1 
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EZjSE 

INSTRTPr=INSTRTPTfl 

GO TO 55 

INDIF 

MSTOT=INSTOITT 



IF (INSTRTPT .B3. 1) THEN 

WRITC (*,*) "Hie first point can only be 1' 

GO TO 215 

WRITE (*,*) 'Uiter first point in waveform file to be processed' 
WRITC (*,*) 'Givai the other parameters chosen thus far,' 

WRTIF (*,*) 'the starting point may range frcni 1' 

WRITE (*,*) ' to',MSTRT 

REM) {*,*) INSTRTPT 

IF (INSTRTPT .GE. 1 .M®. INSTRTRT .LE. MSTRT) TWN 
IF (DSET) GO TO 215 
GO TO 1 
ENDIF 

WRITE (*,*) 'Biter starting point again' 

WRITE (*,*) ' ' 

GO TO 65 
H®IF 

IF (DSET) GO TO 215 

IF {.NOT. DSET) NUFni)=.TRUE. 

IF (.NOT. DSET) NSTRTPT=1 

WRITE {*,*) 'Biter title of data file to be read' 

READ (*,100) iniiE 
0PBN(12,FIIiE=OTLE) 

READ (12, 100) HEADER 
READ (12, 110) NPTS 
IF (NPTS .CT. 1024) THEN 

WRITE (*,*) 'Number of points in data file exceeds the dimension' 

WRITE (*,*) 'of the array used in the program to store the file' 

STOP 

ENDIF 

CL0SE(12) 

IF (NUFILE) THEN 
GO TO 3 

ELSEIF (NSTRTPTf{Kd-m-l)*iaLTAy .LE. NPTS) THEN 

GO TO 215 

ELSE 

GO TO 6 

ENDIF 



IF (NUFILE) THEN 
MAXIMUM=69-Kn-1 



IF (MAXIMUM .GT. NPTS-69) MAXMJHJPTS-69 
MIN=2 

IF (MIN .EQ. MAXIMUM) THEN 
Rd=MIN 

WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE (*,*) 'IGd must be ',MIN 
GO TO 4 
ENDIF 

WRITT (*,*) 'Ihter Kd, >= the estimated order of the system ' 

WRITE (*,*) 'Given the other parameters chosen thus far,' 

75 WRITE (*,*) 'Kd may range from', MIN 

WRITE (*,*) ' to', MAXIMUM 

READ (*,*) Kd 

IF (Rd .GE. MIN .AND. Kd .LE. MAXIMUM) GO TO 4 
GO TO 75 

ELSEIF (DSET) IMEN 
MAXIMUM=41-Kn-1 

IF (MAXIMUM .GT. NPTS-M) MAXIMUM=NPTS-M 

MIN=2 

IWCOaMUM 

85 IF (NSTOTPTf(N+M-l)*DELTAY .LE. NPTS) TOEN 
MAXIMUM^ 

IF (MIN .BQ. MAXIMUM) THEN 

Kd=4tIN 

GO TO 215 

ELSEIF (MAXIMUM .LT. MIN) raEN 
KLTAY=1 

IF (1+(24M-1)*IELTAY .LE. NPTS) TOEN 
Rd=2 

GO TO 135 
ENDIF 

WRITE (*,*) 'Error. Kd must be less than 2' 

Rd=2 

GO TO 215 
ENDIF 

WRITE (*,*) 'Given the other parameters chosen thus far,' 

95 WRITE (*,*) 'Kd may range from ',MIN 

WRITE (*,*) ' to', MAXIMUM 

WRITE (*,*) 'Ehter Kd' 

READ (*,*) Kd 

IF (Kd .(31. MIN .AND. Kd .LE. MAXMM) GO TO 215 
GO TO 95 
EiLSE 
EWhl 
GO TO 85 
ENDIF 
ENDIF 
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C Determine M 

4 IF (NUFILE) TOHJ 

WRTn: (*,*) 'Enter M, the row dimavsic*i of the data matrix' 
IF (.NOT. DSBT .AND. LCNG) 'IffiN 
WRITE (*,*) ' ' 

WRITE (*,*) 'Note: Kd+M points in ', title 
WRITE (*,*) ' will be processed ' 

WRITE (*,*) ' ' 

EJIDIF 

105 WRITE (*,*) 'M may range froD',Kd 
IF (NPTS^M .CT. 69) TTEI 
WRITE (*,*) ' to 69' 

ELSE 

WRITE (*,*) ' to',NPTS-Kd 

EUDIF 

READ (*,*) M 
IF (M .CT. 69) TTEI 

WRITE (*,*) 'M must also be less than 70' 

GO TO 105 

ELSEIF (M .LT. Kd) TTffN 

WRITE (*,*) 'M must be greater than or equal to Kd, Kd= ',Kd 
GO TO 105 

ELSEIF (Kd-W .(?r. NPTS) THEN 

WRITE (*,*) 'Kd+M must be less than or equal to' ,NPTS, ' , ' 
WRITE (*,*) 'the number of data points in',TITTil 
WRITE (*,*) ' ' 

GO TO 105 
EKDIF 

C Begin part for data already set 
ELSE 
N=4Gd 

115 IF (NSTRTPT+(Ka-tN-l)*IffLTAY .LE. NPTS) THEN 
N=N+1 
GO TO 115 
ELSE 
EHM 
ENDIF 

IF (N .BQ. Kd) TEEN 

WRITE (*,*) 'M must equal' ,Kd 

EH(d 

GO TO 215 

EUDIF 

MAXIMUM=N 

IF (MAXMJM .GT. 69) MAXIMUW=69 
IF (Kd+Kn+1 .E3. MAXIMUM) TEEN 
M=Kd+Kn+l 
GO TO 215 

ELSEIF (Kd+Kn+1 .CT. MAXIMUM) TEEN 
WRITE (*,*) 'Kd must be reduced' 
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GO TO 3 
ELSE 

MDHCd-HCn+1 

INDIF 

IF (MIN .LT. Kn+Kd+1) MDMn+Kd+1 
125 WRITE (*,*) 'M may range fran',MIN 
WRITE (*,*) ' to’, MAXIMUM 

WRITE (*,*) ’alter M’ 

READ (*,*) M 

IF (M .GE. MIN .AND. M .LE. MAXIMUM) GO TO 215 

GO TO 125 

ENDIF 



c 

135 

5 

145 
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Determine DELTAY 

IF (.NOT. NUFILE) GO TO 215 

N=1 

IF (NSTRTET+N*(Kd4M-l) .LE. NPTS) TEIN 
N=NH 
GO TO 145 
ELSE 

ENDIF 

IF (N .BQ. 1) TEEN 

WRITE (*,*) 'Given the other parameters chosen thus far,’ 
WRITE (*,*) ’Spacing can oily be 1’ 

EEiTAY=l 

IF (NUFILE) THEN 
GO TO 165 
ELSE 

GO TO 215 
ENDIF 
ENDIF 



IF (.NOT. DSET .AND. LCNG) TEEN 
WRITE (*,*) ’Eiiter spacing between the ’,Kd+M 
’data points of ’, TITLE 
’to be processed ’ 



WRITE (*,*) 
WRITE (*,*) 



WRITE (*, 

(*, 



*) 

*) 

*) 

*) 



WRITE 
WRITE (* 

WRITE (* 

WRITE (*,*) 

ELSE 

WRITE (*,*) 

WRITE (*,*) 

ENDIF 

WRITE (*,*) 

WRITE (*,*) ’ 

READ (*,*) MliTAY 
IF (DELTAY .GE. 1 .AND. 
IF (NUFILE) TEEN 



’If, for example, c»ie is chosen, thoi ’,Kd+M 
’consecutive points in ’,’nTLE 
’will be processed ’ 



’Biter spacing 



’Spacing may range fran 
’ to’,N 



MLTAY .LE. N) THEN 
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GO TO 165 

GO TO 215 
IMIIF 
W.SF: 

GO TO 155 
H®IF 

165 WRITE {*,*) 'Do you wish to adjust eigenvalues? (y/n) ’ 

READ (*',150) WI 

IF (YN .IQ. 'N' .OR. YN .BQ. 'n') IBEN 
IF (NUFILE) GO TO 6 
GO TO 215 
ENDIF 

IF (YN .NE. 'Y' .AND. YN .NE. *y') GO TO 165 
2 WRITE (*,*) 'Discard or canpensate eigenvalues? (d/c) ' 

READ (*,150) DC 

IF (DC .EQ. 'D' .CR. DC .B2. 'd') TBOI 

NRT=W 

GO TO 175 

ENDIF 

IF (DC .NE. 'C .AND. DC .NE. 'c') GO TO 2 

WRITE (*,*) 'Biter estimate of the actual order of the system' 

WRITE (*,*) ' ' 

IF (LGM3) THEN 

WRITE (*,*) 'This estimate will be used to determine the ' 
WRITE (*,*) 'number of eigenvalues compensated or discarded ' 
ENDIF 

175 WRITE (*,*) 'the estimate may range frcm 2' 

WRITE (*,*) ' to',Kd+Kn+l 

READ (*,*) NRT 

IF (NRT .GT. Kd+Kn+1 .CR. NRT .LT. 2) THEN 
GO TO 175 

EliSEIF (.NOT. NUFILE) THEN 

GO TO 215 

ENDIF 

6 NSTRTPT=1 

185 IF (NSTRTPT+(Kd4M-l)*IMiTAY .LE. NPTS) THEN 
NSTTm*WISTHTPT4-l 
GO TO 185 
ELSE 

NSTRTPW(STHTPT-1 

ENDIF 

IF (NSTRTTT .EQ. 1) THEN 

WRITE (*,*) '(Jiven the other parameters chosen thus far,' 

WRITE (*,*) 'the starting point for processing the data' 

WRITE (*,*) 'must be the first point in the data file' 

GO TO 215 
ENDIF 
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WRITE (*,*) 'filter desired starting point in data file' 

IF (.NOr. DSET .AND. LONG) THEN 

WRITE (*,*) '1 indicates the first point in the data file ' 
BNDIF 

WRITE (*,*) ' ' 

WRITE (*,*) 'Given the other parameters diosen thus far, ' 
195 WRITE (*,*) 'the starting point may range fran 1' 

WRITE (*,*) ' to',NSnOTT 

READ (*,*) N 

IF (N .GE. 1 .AND. N .LE. NSTRITT) THIN 

NSTRTET=N 

ELSE 

WRITE (*,*) 'filter starting point again' 

WRITE (*,*) ' ' 

GO TO 195 
INDIF 

IF (.NOT. NUFILE) GO TO 215 
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IF (DSET) THIN 

IF (NCAUS .EC. 1) THEN 

NCAUS=2 

GO TO 215 

ELSE 

NCAUS=1 

GO TO 215 

INDIF 

INDIF 

WRITE (*,*) 'Do you want the data matrix arrangement to tie' 
WRITE (*,*) ' ' 

WRITE (*,*) '1. Causal' 

WRITE {*,*) '2. Non-causal' 

WRITE (*,*) ' ' 

WRITE (*,*) 'Please alter 1 or 2 ' 

READ (*,*) NCAUS 
IF (NCAUS .EC. 1) THEN 
CAUSAL=.THUE. 

ELSEIF (NCAUS .EC. 2) THEN 
CAUSAL=.FALSE. 

GO TO 205 
INDIF 
GO TO 215 



12 WRITE (*,*) 'Enter title of file to contain parameters' 
READ (*,100) TTIL 
CPEN(l,FIIi>=nTL) 

WRnE(l,100) TTILE 
WRnE(l,110) NPTS 
WRITE (1,110) NRT 
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WRITE(1,U0) Kd 

wm: (1,110) M 

TOTTC (1,110) EELTAY 

vRrn:(i,iio) nstoitt 

WRITE (1,110) NCAUS 
WRITE(1,100) TTILD 
WRTn: (1,110) NDPTS 
WRITC (1,110) Kn 
WRITE(1,110) DISTRTOT 
CLOSE(l) 

IF (DSET) GO TO 215 

15 IF (DSET) TOEN 
ClOSE(2) 

CLOSEO) 

CALL SUBPLT(NCfVIFLAY) 
B©IF 



215 DSET=.TRUE. 

wRTn:(*,*) ' ' 

WRrn:(*,*) ’l. Data file to be processed ',T 

+ITTJE 

WRTn:(*,*) ’ Number of data points in data file ',NPTS 

WRrra;(*,’*) ’2. Estimated order of the system ',NRT 

WRTn;(*,*) ’3. Kd, the number of columns in the data matrix', Kd 
WRTITK*,*) *4. M, the number of rows in the data matrix’, M 
WRITE(*,*) '5. Spacing between data points being processed ',IMjTA 
+Y 

WRTITK*,*) *6. First point in the data file to be processed' ,NSTRT 
+PT 

WRTITK*,*) ’ Last point in the data file to be processed’ ,NSTRT 
+PT+Kd4M-l 

IF (NCAUS ,BQ. 1) 'HEJ 

WRrn:(*,*) Data matrix arrangemait for processing CA 

-KJSAL ’ 

ELSE 

WRITE (*,*) ’7. Data matrix arrangemait for processing M3N-CA 

■HJSAL ’ 

H©IF 

WRITE(*,*) ’ ’ 



WRTIE(*,*) ’8. File containing excitation waveform ’,T 

+nijD 

WRTIE(*,*) ’ Number of data points in above file ', NDPTS 

WRTIE(*,*) ’9. Estimated order of the waveform ',Kn 

WRTIE(*,*) ’10, First point in the file to be ’ 

WRTIE(*,*) ’ input into the data matrix ’,INSTR 

+TPT 
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VRITE(*,*) 



t t 
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VRITE(*,*) 'll. Begin processing using above settings' 

VRITE{*,*) '12. Store parameters 1-10 in a file' 

WRITE (*,*) '13. Retrieve parameters 1-10 from a previously created 
+ file' 



WRrTE(*,*) '14. Reset overlays' 

WRrrE(*,*) '15. Re-plot overlays' 

WRITE{*,*) '16. Bid this session of Cadawf-Solcccn signal process! 
■mg' 

WR1TE(*,*) ' ' 

WRITE (*,*) 'Enter an integer frcro 1 to 16 to make changes as oftai 
+ as you desire' 

READ (*,*) NMENU 

IF (NMEMJ .LT. 1 .OR. MENU .CT. 16) 'lEEN 
WRTIE(*,*) 'Biter an integer from 1 to 16' 

GO TO 225 
HTOIF 



GO TO (1,2, 3, 4, 5, 6,7, 8, 9,10,11, 12, 13, 14,15, 16), NMENU 

11 CPEN(12,FILE==nTLE) 

RERD(12,100) HEMM 
REaD(12,110) NPTS 
RERD(12,120) XQ 
READ(12,120) XQ 
DO 235 I=1,NPTS 
RERD(12,120) D(I) 

235 CONTINUE 
CLOSE (12) 

C5TN (8 , FILE=TITLiD) 

RERD(8,100) HEADER 
READ (8, 110) NDPTS 
READ (8, 120) 

READ(8,120) 

DO 245 I=l,NIFrS 
READ(8,120) Dx{I) 

245 CCNTINUE 
CL0SE(8) 

KdPLT:^ 

WRriE(*,*) 'enter title of file to contain real part of poles' 
READ(*,100) TriLER 
OPEN (2 , FILE)=nTIjER) 



WRITE {*,*) 'enter title of file to contain imaginary part of poles' 

READ(*,100) TITIiEI 

CPEN(3,FTLE=nTLEI) 
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WRITE (10, 130) (KdPLT) 
WRITE (10, 1(X)) TnUR 
WRITE (10, 100) TnUEI 
130 FORMAT (12) 

>W=MAX(M,Kd+Kn+l) 

100 FORMAT (A) 

no FORMAT (15) 

120 FORMAT (E12. 6) 

150 FORMAT (A) 



DO 255 I=l,Kd4M 

Dy(I)=®( (I-l)*EELTAY+iKITRTPT) 

255 (XOTINUE 

265 DO 285 1=1, M 

DO 275 J=l,Kd+Kn+l 
A(I,J)=Dy(I+J) 

IF (J .GE. Kd+1) A(I,J)^(I4J+INSTRTPT-2-Kd) 

275 OOWTINUE 
285 OCNITNUE 

B(l)=Dy(l) 

DO 295 1=2, M 
B(I)=A(I-1,1) 

295 CONTINUE 

N^-Hto+l 

C Begin singular value decomposition 

CALL SVD(MACHEP,M,N,MN,A,W,MA1U,U,MATV,V,IIRR,RV1) 

C Errors in SVD? 

IF dm .CT. 0.0) TEEN 

WRITE (*,*) 'Error in singular value number ',im,ST:P 
INDIF 

IF (YN .EQ. 'N') GO TO 385 

DO 305 I=l,Kd4Kn+l 
XP(I)=0.0 
305 CONTINUE 

C Discard or conpoisate eigenvalues 

c Order singular values 

XP (1)^(1) 
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DO 335 I=2,Kd+Kn+l 
DO 325 C^1,I 

if (W(I) .GT. XP(J)) 'IBEN 
DO 315 K=I+1,J,-1 
315 XP(K)=XP(K-1) 

XP(j)=W(i) 

GO TO 335 
I«DIF 

325 ocrnNUE 
XP(I+1)=^(I) 

335 CXWnWUE 

C XP( ) now contains ordered singular values: XP(1) is the largest 

C Discard eigenvalues 

IF (DC .EQ. 'D') TOEN 
DO 345 J=4®T+l,Kd-HCn+l 
345 W(J)=(0.0) 

ELSE 

C Cciiipensate eigenvalues 

AVG=0.0 

DO 355 J=NRTM-l,Kd-HCn+l 
AVG=AVO+XP(J)**2 
355 CCNimJE 

IF (Kd+Kn+1 .GT. NRT) AV(>=AVG/DBIi:(Fl;OAT(Kd-HCn+l-NRT) ) 

DO 375 J=l, Kd+Kn+1 
DO 365 K=l, Kd+Kn+1 
IF ( W(J) .EQ. XP(K) ) TOEW 
IF ( K .GT. NRT ) THEN 
W(J)=0.0 
ELSE 

H ( J) =DS(3^ (DABS ( W ( J) *W ( J) -AVG) ) 

ENDIF 
GO TO 375 
ENDIF 
365 OCNTINUE 
375 OCNTINUE 
ENDIF 



385 DO 405 1=1, M 
DO 395 J=1,M 
irr(I,J)=(U(J,I)) 

395 CONTINUE 
405 CONTINUE 

C Form SK3MA+ (Kd+Kn+1 x M) 
DO 425 1=1, Kd+Kn+1 
DO 415 J=1,M 
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SI(31R(I,J)=0.0 

IF (I .EQ. J .AND. W(J) .NE. 0.0) TOEN 
SIGMA(I,J)=1.0d0/W(J) 

SIGMA(I,J)=O.ODO 

ENDIF 

415 OCKTINUE 
425 OCWniWE 

C Form SIGMA (M x Kd-HCn+1) 

DO 445 1=1, M 
DO 435 J=l,Kd+Kn+l 
SIG(I,J)=0.0 

IF (I .EQ. J) SIG(I,J)=W(J) 

435 OCMTINUE 
445 OCNTINUE 

C V=4(d+Kn+lxRa4fti+l,SIGMA-H=i03-HO^ 

CALL MXMUL(V,SIGMA,Kd+KiH-l,Kd-H(n+l,M,VS) 

C VS=Kd4Kn+lxM,irr=Mrfl,AINV^ 

CALL MXMUL(VS,Ur,Kd+Ktt+l,M,M,AINV) 

C Calculate matrix multiplication of AINV x B, where 

C AIN\^4fOYflrfl,B=4ixl,XP=«d4Kn^ 

CALL MXMUL(AINV,B,Kd-H(n+l,M,L,XP) 

C Canpute autoregressive coefficients from prediction coefficients 

IF (XP(Kd) .EQ. 0.0) THEN 

write; (*,*) 'E3UIOR, avoiding division by zero' 

STCP 

ELSE 

B(Kd)=1.0d0/XP(Kd) 

EKDIF 

DO 455 1=2, Kd 
B <1-1 ) :^B (Kd) *XP (Kd-i+1 ) 

455 OCWnNUE 

DO 465 i=l,Kd 
X(I)=-B(Kd-I+l) 

IF (NCAUS .EQ. 1) X(I):^XP(Kd-I+l) 

465 CCOTINUE 

X(Kd+l)=1.0 

C Conpite the roots of the polyncndal in z 

CALL PCm’(X,CXF,KD,ROCriR,ROCm,IER) 

IF (m .NE. 0) WRITE (*,*) ’E3(R(* with PCMT, IEK=' ,IIR,STCy 
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DO 475 1=1 ,Kd 
WRITE(2,120) ROOTO(I) 

WRrrc(3,120) ROOTI(I) 

S (I) =<)CMPIiX(ROOTR (I) , I«DOTI (I) ) 

475 CCNmWE 

MA(3<Xi=0 
DO 485 1=1, Kd 

IF (CDABS(Sd)) .GE. l.ODO) MAGPC*>=tlAGPa/<-l 
485 OCNTINUE 

VRITE(*,*) 'f of poles vith magnitude <= 1’ ,Kd-MAGP(Xi 
WRITE (*,*) 'HIT ANY KEY TO CXmUUE’ 

READ (*,100) HEAIER 

C Plot poles 

»7VH^Y=NOVimY+l 

CL0SE(2) 

CLOSED) 

CALL SUBPLT(NOVIPLAY) 

J=0 

K=0 

DO 495 1=1, Kd 

IF (CDABS(S(D) .LT. 1.0) THEN 
WRITE (*,*) S(I),CDABS(S(I)) 

J=J+1 

K=K+1 

ENDIF 

IF (J .EQ. 20) TOEN 

WRITE (*,*) 'HIT ANY KEY TO CONTINUE' 

READ (*,100) HEAMR 

J=0 

ENDIF 

495 OCNTINUE 

WRITE (*,*) 'Poles with magnitude less than one ',K 
WRITE (*,*) 'HIT ANY KEY TO CONTINUE' 

READ (*,100) HEftDEN 

GO TO 215 

16 STCy 
El® 

*Z 
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APPENDIX C. MATRIX MULTIPLICATION 



SUBROUTINE MXMUL(A,B,RA,CA,CB,AB) 

DTIBGro RA,CA,CB 

REAL*8 A(70,70),B(70,70),AB(70,70) 

C Calculates matrix multiplication of A x B^, i4iere 

C A=RAxCA,B=CAxCB,ABi:^^AxCB 

DO 30 1=1, RA 
DO 20 J=1,CB 
AB(I,J)=0.0 
DO 10 K=1,CA 

AB{I,J)=AB(I,J)+A(I,K)*B(K,J) 

10 OCMITNUE 

20 CCKTINUE 

30 CONTINUE 

RETURN 

END 
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APPENDIX D. GRAPHICS ROUTINE 



SUBROUITNE SUBPLICNOVERIAY) 



C 

C MS-FC*TOAN Program using "Grafmatic" Library Subroutines. 

C Plots a Solid Line and Optional Overlay Plot for Ccopariscn. 

C Writtai ty M.A. Morgan with Latest Update August 1989. 

C 

C Default Printer is "IBM Graphics" (e.g. E^son, (Mdata, IBM) 

C With Plot Rotated 90 degrees Fran the Vertical. "GrafPlus.Com" 

C May be Run to Rotate Plot Upright on Paper and to Use a Variety 

C of Intact Printers. "<a-af Laser. Com" May be Run to Use a Laser 
C Printer. See (hraf Plus/Laser Manual Frcm Jewell Ttechnology. 

C 

C 

CHARACUR*! YN,YNl,DUM,YN2,SYMB(X,BELL,FEED,FFyN 
CHARACITR*4 LINE 
CHARACTER*? SYMB 

CHARACTER*16 LTCT,CTn,ENAME,TTniR,TriliE3 
CHARACTER*64 TTHiE.HOCfY 
REAL CRTR(70),CRTK70),NRTR(70),NRn(70) 

INTE)t3R*2 N,JRCW,JCCi,ISYMl,ISYM2,mPEl,ITYPE2,NSCRN 
INTE)GER*2 CYAN,(»EIN,WHITE,YELLCW,RED, BLACK, BLUE, EnW) 

IWrEX3R*2 JRC«l,JRC«2,J0CXil,JCC4i2, CROSS, KdPLT,I 
INTBGER*2 PURPLE,RUST 
EmRNAL XEUN,YFUNP,YFWN 
UNE=' — ' 

WHriT)=7 

GREm=10 

CYAN=11 

YELI£W=14 

RED=12 

BLACK=0 

BLUE=1 

Enw>2 

PURPLE>=5 

RUST=6 

BELL=CHAR(7) 

FEED<HAR{12) 

C Clear Screen and Put Introduction - on Blue Backgcwnd for E)GA 
C Only; Another Bac)cground Color is Possible by Changing "BLUE" 

C in the Calls to and QOVSCN. 

CALL QSM»E(NTWD) 

CALL ^REG(0,BLUE) 
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CWIi QCVSCN(BLUE) 

WRITE(*,*) BELL 
NS=1 

NSCRN =16 
mPE2=0 

C Calling GRAFMATIC Routines and Plotting FI Solid Line Graph 
ITYPEIO 
ISYM1=-1 
NDorsi=o 

JRC«1=1 
JRCW2=350 
JCm= 75 
J00L2= 565 
XMIN=-1.2 
XMRX=1.2 
YMDt=-1.2 
YMaX=1.20 
TOVm=1.115 
XC*G=0.0 
YDRG=0.0 
XST^l.l 

xnN=i.i 

YST^l.l 

YFIN=1.1 

25 CALL QSMCOS (NSCRN) 

CALL (^>LCrr(JCCU,JC(^,JR(m,JROW2,XMIN,XMAX,Y>aN,Ym,}OM(G,ra^ 
+1,Y0VERX,1.5) 

CALL QSEnUP(NDOTSl, CYAN, ISYMl, RED) 

IF(XFIN-XST .LE. 9.0) XMAJ0R=0.6 
IF(XFIN-}CT .LE. 6.0) XMAJ(»=0.4 
IF(XnN-XST .LE. 3.3) XMAJC»=0.2 
IF(XFIN-XST .GE. 9.0) XMAJCR=(XFIN-XST)/10.0 

MINC*=0 

LABEL=1 

NDBC=2 

CALL Q50DaS(XST,XFIN,XMAJ(»,MIN0R,LABEL,NKi:) 

YMAJC«=XHAJCR 

CALL gn«IS(YST,YFIN,YMAJ(»,MINC», LABEL, NDTC^ 
c Plot unit circle 
A=-1.0 
B=1.0 

CALL QOJRV(XFUN,YFUNP,A,B) 

CALL QCURV(XFWl,YFtJNN,A,B) 

IF (NOViRLAY-1 .LT. 1) THEN 
IF (NCVERLAY-1 .EQ. 0) THEN 
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WRITE (*,3) NDVERIAY-1 

ELSE 

NZE3^^^ 

WRITE (*,3) NZERO 
ENDH 

ELSEIF (NCWMAY-l .CT. 1) THS^ 

WRITE (*,3) NOVERLAY-1 

ELSE 

WRITE (*,4) NDVERIAY-1 
E1®IF 

3 FORMAT {13/ OVERLAYS') 

4 EDRMAT (13/ OVERLAY ’) 

REWIND (10) 



DO 20 I=l,NOVERLAY 
READ (10,110) KdPLT 
READ (10,100) TITLER 
READ (10,100) TITLEI 
aPEN{2,ETLE=nTLER) 

CfOI{3,FILEMTITiEI) 

NKd=KdPLT 
DO 27 J=l, KdPLT 
READ (2,120) NRTR(J) 

READ (3,120) NRTI(J) 

IF (DSQRT{NRTO{J)**2-tNRTI{J)**2) .CT. 1.1) THEN 

NKd=NKd-l 

NRTR(J)=0.0 

NRn(J)=0.0 

ENDIF 

27 CCNTINUE 
FURPLB=5 
RUST=6 
WHITE=7 
(REEN=10 
CYAI^ll 
YELLCW=14 
RED=12 
BLUE=1 

IF (I .EQ. 1) 'ffiEN 
CALL QSETUP{ND0TS1,CYAN,ISYM1,RED) 

ELSEIF (I .EQ. 2) THEN 
CALL QSETUP(NDOrSl,CYAN,ISYMl,GREEK) 

ELSEIF (I .EQ. 3) TOER 
(ALL QSEnUP(NDCrrSl, CYAN, ISYMl, YELLOW) 

ELSEIF (I .EQ. 4) 71®! 

CALL QSEnJP(NDOTSl, CYAN, ISYMl, BLUE) 

ELSEIF (I .EQ. 5) niER 
CALL QSEIUP(NDCrrSl,CYAN,ISYMl,WHITE) 

ELSEIF (I .EQ. 6) THER 
CALL QSEIUP(NDOrSl, CYAN, ISYMl, PURPLE) 
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20 

100 

110 

120 

160 

40 



HLSEIF (I .BQ. 7) 'IHEU 
CALL QSEIUP (NDCrrSl, CYAN, ISYMl, RUST) 

ELSE 

CALL QSEnJP{NDOTSl,CYAN,ISYMl,RiI)) 

ENDIF 

CALL giML{rryPEl,KdPLT,NRTO,NRn) 

OCNTINUE 

READ (*,100) DUM 
GO TO 40 

HCOPY=’ HARDCOPY — > HTTER P p* 

CAUi ^yTJrr(30,H0CPy,RED,25,l) 

CAIi QCM0V(55,1) 

HCCPY=’ 

CAIi QpTCT(40,H0CPY,BLACK,25,l) 

IF (DUM .NE. 'P' .AND. DUM .NE. ’p’) GO TO 40 

CALL CPSCRN 

OPEN (1,FILE='PRN’) 

WRITE (1,160) FEED 
FORMAT (A) 

F0RMAT(I2) 

FO»IAT(E12.6) 

FORMATC ’,A,\) 

CONTINUE 

CALL QSEEl®(NTWD) 

CALL (yREG(0,BLUE) 

CALL QCVSCN(BLUE) 

WRriE(*,*) NKd, 'points were plotted' 

RETURN 

END 

REAL EimiON XFUN(T) 

XEUN=T 

RETURN 

END 

REAL FUNCTION YFTM»(T) 

YEUNP=S(^T(1.0-T*T) 

RETURN 

END 

REAL FUNCTION YFUN(T) 

YFUNN=-Sg?r{1.0-T*T) 

RETURN 

END 
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